Naval  Research  Laboratory 

Washington,  DC  20375-5320 


NRL/MR/8233— 05-8897 


Investigations  of  Spacecraft  Orbits 
Around  the  L2  Sun-Earth  Libration  Point 

Alan  M.  Segerman 

AT&T  Government  Solutions,  Inc. 

Vienna,  VA 


Michael  F.  Zedd 

Control  Systems  Branch 
Spacecraft  Engineering  Department 


July  29,  2005 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway, 
Suite  1 204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 
29-07-2005 

2.  REPORT  TYPE 

Memorandum  Report 

3.  DATES  COVERED  (From  -  To) 

January  2002-September  2004 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

NNG04CB011 

5b.  GRANT  NUMBER 


Investigations  of  Spacecraft  Orbits  Around  the  L2  Sun-Earth  Libration  Point 

5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION  REPORT 

NUMBER 


6.  AUTHOR(S) 

Alan  M.  Segerman*  and  Michael  F.  Zedd 


Naval  Research  Laboratory,  Code  8233 
4555  Overlook  Avenue,  SW 


Washington,  DC  20375-5320 


NRL/MR/8233— 05-8897 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR  /  MONITOR’S  ACRONYM(S) 


NASA/GSFC 


Goddard  Space  Flight  Center 
Greenbelt,  MD  20771 


11.  SPONSOR /MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 

*AT&T  Government  Solutions,  Inc.,  Vienna,  VA 

This  Memorandum  Report  is  a  compilation  of  three  separate  reports,  of  increasing  complexity,  for  the  sponsor  (none  previously  published). 

14.  ABSTRACT 

Space  agencies  are  planning  missions  to  the  vicinity  of  the  Sun-Earth  L2  point,  some  involving  a  distributed  system  of  telescope 
spacecraft,  configured  in  a  plane  about  a  huh.  An  improved  understanding  is  developed  of  their  relative  motion.  First,  the  telescope 
equations  of  motion  are  written  relative  to  L9  in  the  context  of  the  classical  circular  restricted  three-body  problem,  and  expanded  in  terms 
of  the  distance  from  L9.  A  basic  examination  is  presented  of  the  spacecraft  configuration  requirements  and  effects  of  small  orbit  insertion 
errors.  Next,  the  telescope  equations  of  motion  relative  to  the  hub  are  written  and  further  expanded  in  terms  of  the  hub-L9  and  hub-telescope 
distances.  An  analytical  solution  is  developed  and  a  halo  telescope  orbit  investigated,  with  appropriate  initial  conditions.  Then,  the  force 
model  is  extended  to  include  perturbations  to  an  accuracy  of  10  to  20  m,  cast  as  additive  contributions  to  the  circular  restricted  problem. 
Perturbations  include  Earth’s  orbital  eccentricity,  lunar  motion,  solar  radiation  pressure,  and  small  thrusting  forces.  Simulations  are 
presented,  along  with  solution  sensitivity  to  errors  in  hub  position  knowledge. 


15.  SUBJECT  TERMS 


Sun-Earth  libration  point;  L2  point;  Constellation;  Formation  flying;  Relative  motion;  Chief;  Deputy;  Hub;  Telescope 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

OF  ABSTRACT 

OF  PAGES 

Michael  F.  Zedd 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

UL 

178 

19b.  TELEPHONE  NUMBER  (include  area 

Unclassified 

Unclassified 

Unclassified 

code) 

(202)  404-8337 

1 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Contents 


Part  1:  Introduction  to  the  Relative  Motion  of  Spacecraft  About 
the  Sun-Earth  X2  Point  1 

Part  2:  Linear  and  Quadratic  Modelling  and  Solution  of  the 
Relative  Motion  43 

Part  3:  Modelling  the  Perturbations  —  Elliptical  Earth  Orbit, 
Lunar  Gravity,  Solar  Radiation  Pressure,  Thrusters  97 


iv 


PREFACE 


This  Memorandum  Report  consists  of  a  compilation  of  three  individual 
reports,  of  increasing  complexity,  describing  investigations  of  formation  flight 
of  spacecraft  in  the  vicinity  of  the  L2  Sun-Earth  libration  point.  The  indi¬ 
vidual  reports  form  the  following  parts  of  this  compilation: 

•  Introduction  to  the  relative  motion  of  spacecraft  about  the  Sun-Earth 
L2  Point 

•  Linear  and  quadratic  modelling  and  solution  of  the  relative  motion 

•  Modelling  the  Perturbations  —  Elliptical  Earth  Orbit,  Lunar  Gravity, 
Solar  Radiation  Pressure,  Thrusters 

The  three  parts  are  self-contained,  with  somewhat  varying  notation  and 
terminology. 

This  work  was  funded  by  the  Distributed  Spacecraft  Technology  Program 
at  NASA’s  Goddard  Space  Flight  Center.  The  sponsor  was  Dr.  Jesse  Leitner 
of  the  Mission  Engineering  and  Systems  Analysis  Division. 

After  fairly  significant  literature  searches,  this  new  work  (of  Parts  2  and  3) 
is  deemed  to  be  unique  because  it  describes  the  primary  perturbations  to  the 
description  of  relative  motion  between  nearby  spacecraft.  The  effect  of  the 
elliptical  motion  of  the  Earth  about  the  Sun  was  verified  to  be  the  dominant 
perturbation  to  the  circular  restricted  three  body  problem.  Contributions 
due  to  lunar  gravity  and  solar  radiation  pressure  are  seen  to  have  much 
smaller  effect. 


v 
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1  Introduction 


The  long-term  goals  of  our  sponsor  are  to  develop  high-fidelity  equations 
of  motion  representing  orbiting  spacecraft  near  the  Sun-Earth  L2  point  and 
equations  of  relative  motion  between  orbiting  spacecraft  near  the  L2  point. 
These  equations  will  be  used  to  develop  orbit  control  schemes  for  a  constel¬ 
lation  of  spacecraft  near  this  L2  point. 

This  report  provides  a  starting  point  to  the  future  development  of  high- 
fidelity  equations  by  first  developing  equations  of  motion  for  the  circular 
restricted  three-body  problem.  Additionally,  equations  are  derived  for  ana¬ 
lytical  expressions  for  the  relative  motion  between  spacecraft  orbiting  near 
the  L2  point. 

We  follow  the  terminology  of  Hamilton’s  thesis  [1]  as  we  present  four 
sets  of  equations  of  motion:  three  nonlinear  sets  and  a  linearized  set.  Their 
derivations  are  in  separate  appendices.  These  equations  were  developed  to 
examine  relative  errors  of  the  linearized  and  quadratic  effect  models  with 
respect  to  the  baseline  set  full  nonlinear  equations. 

The  four  types  of  models  for  the  circular  restricted  three-body  problem: 

•  full  nonlinear  model  of  spacecraft  motion  with  respect  to  the  system 
barycenter 

•  full  nonlinear  model  of  spacecraft  motion  with  respect  to  the  L2  point 

•  linearized  model  of  spacecraft  motion  with  respect  to  the  L2  point 

•  quadratic  (nonlinear)  model  spacecraft  motion  with  respect  to  the  L2 
point. 

We  formed  a  Taylor  series  expansion  of  the  nonlinear  equations  about 
the  libration  point.  Our  linearized  and  quadratic  models  are  respectively  the 
first-  and  second-order  Taylor  series. 

The  analytical  expressions  of  relative  motion  are  based  on  the  solutions 
to  the  linearized  equations. 

2  Problem  Definition 

The  restricted  three-body  problem  is  a  simplification  to  the  general  three- 
body  problem.  The  general  three-body  problem  is  the  description  of  orbital 
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motion  of  three  bodies,  of  arbitrary  mass,  in  arbitrary  orbits.  The  general 
three-body  problem  has  not  been  solved  in  closed  form.  “Work  has  focused 
on  simplifying  the  general  problem.  One  special  analytical  solution — the 
restricted  three-body  problem — has  been  known  since  the  time  of  Euler  and 
Lagrange.”  [2] 

“We  define  our  problem  as  follows:  Two  bodies  revolve  around  their  cen¬ 
ter  of  mass  in  circular  orbits  under  the  influence  of  their  mutual  gravitational 
attraction  and  a  third  body  (attracted  by  the  previous  two  but  not  influenc¬ 
ing  their  motion)  moves  in  the  plane  defined  by  the  two  revolving  bodies. 
The  restricted  problem  of  three  bodies  is  to  describe  the  motion  of  this  third 
body .”  [3] 

The  two  revolving  bodies  are  called  the  primaries.  The  masses  rri \  and 
m2  of  these  bodies  are  arbitrary  and  are  considered  as  point  masses.  The 
mass  of  the  third  body,  m3,  is  infinitesimal  and  does  not  influence  the  motion 
of  masses  mi  and  m2.  In  our  problem,  consider  the  m3  mass  to  be  a  space¬ 
craft.  Also  consider  the  sum  of  the  mass  of  the  Earth  and  its  Moon  to  be 
m2.  Assume  the  Earth  and  Moon  are  in  circular  orbit  about  their  common 
barycenter  (system’s  center  of  mass).  In  turn,  this  barycenter  is  in  circular 
orbit  about  the  Sun,  which  is  mass  m  1 .  These  four  masses  compose  our 
system;  however,  there  are  two  large  bodies.  The  larger  of  the  large  bodies 
is  the  Sun  and  the  smaller  of  the  large  bodies  located  at  the  Earth-Moon 
barycenter  (mi  >  m2). 

We  copy  portions  of  Hamilton’s  thesis  Section  2.1  to  define  terminology 
and  coordinate  systems. 

For  this  preliminary  analysis,  the  orbits  of  the  Moon  around  the  Earth 
and  the  Earth  around  the  Sun  are  assumed  to  be  circular,  have  a  constant 
angular  speed,  and  be  in  the  same  orbital  plane.  (Spacecraft  motions  are  not 
restricted  to  the  plane.)  When  these  assumptions  are  made,  the  restricted 
three-body  problem  becomes  the  more  specific  circular  restricted  three-body 
problem. 

The  rotating  coordinate  frame  shown  in  Figure  1  is  used  for  the  devel¬ 
opment  of  the  equations  of  motion.  The  origin  is  the  system  barycenter  of 
mi  and  m2.  The  first  basis  vector,  d\,  points  toward  the  smaller  body.  The 
third  basis  vector,  03,  is  parallel  to  the  orbit  normal  with  origin  at  the  system 
barycenter.  The  second  basis  vector,  a2,  is  defined  by  the  cross  product  of  a3 
and  ai  (0,2  =  (13  x  ai).  These  vectors  form  an  orthogonal  coordinate  frame. 

The  position  vector  of  the  spacecraft  in  this  rotating  frame  is 


6 
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Figure  1:  The  Three-Body  Problem’s  Rotating  Coordinate  System. 


r  =  Xai  +  Ya  2  +  Za 3. 

Writing  vectors  from  the  large  bodies  to  spacecraft  in  the  rotating  frame 
gives 


f\  —  ( X  +  D\)ai  +  Y  a2  +  Zeis 

and 

f2  =  (A  —  D2)ai  +  U02  +  Zeis, 

where  D 1  is  the  distance  from  the  barycenter  to  rn  1  and  D2  is  the  distance 
from  the  barycenter  to  m2-  Note  that 


D  —  D\  +  D2. 


Use  the  equation  for  the  definition  of  center  of  mass:  rn  \  D 1 
calculate  the  magnitudes  of  l)  1  and  D2  given  D: 


Dl 

D2 


D 


1  +  m 

m2 

D 


1  + 


TTL2  ’ 
mp 


m2D2  to 
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Define  the  constant  angular  speed  of  the  frame  about  03  as 


lG(m  1  +  m2)  _  I G(m sun  +  mEarth  +  rnMoon ) 

U  ~  V  D3  “  v  D 3 

where  G  is  the  gravitational  constant.  As  shown  in  Appendix  A,  we 
define  the  constant  angular  speed  of  the  frame  about  03  as  the  two-body 
mean  motion  Cj  =  ua^ 

Also  define 


l-i  t  =  Grrii  =  GmSun 


and 

H 2  G / / / 2  G(m Earth  A  UCv/oon) * 

With  these  definitions,  we  state  four  sets  of  equations  that  describe  the 
motion  of  the  spacecraft.  The  derivation  of  these  equations  is  shown  in 
respective  appendices. 


3  Full  Nonlinear  Equations  for  Circular  Re¬ 
stricted  Three-Body  Problem  (relative  to 
barycenter) 

These  equations  describe  the  spacecraft’s  motion  relative  to  the  system’s 
barycenter  in  a  rotating  coordinate  frame.  The  frame  rotates  with  a  constant 
angular  speed. 


X  -2uY  -uj2X  = 
Y  +  2uX-uj2Y  = 
Z  = 


/Mi(X  +  D\ )  /i2(  X  —  D2 ) 


N3 

ViY  H2Y 


1^2 1 ; 


n 


r  2 


Hi  Z 


|d|: 


'^2 


(1) 

(2) 

(3) 


where  |r*i  |  and  |r2|  are  the  magnitudes  of  the  vectors  rq  and  r2,  respectively. 


(not  to  scale) 


Figure  2:  Libration  Point  Locations. 


The  initial  conditions  to  be  specified  are 

A'(0)  A'(0) 

m  y( o) 
m  z( o). 

These  equations  are  derived  in  Appendix  A. 

By  setting  all  time  derivatives  in  these  equations  of  motion  to  zero,  five 
libration  points  can  be  calculated.  The  location  of  these  points  depend  on  the 
masses  and  distances  of  the  bodies;  however,  three  points,  shown  in  Figure  2, 
are  always  collinear  with  the  two  large  bodies  (Li,  L-2,  L3),  and  the  other 
two  points  form  equilateral  triangles  with  the  two  large  bodies  (L4  and  L5). 

The  motion  near  the  collinear  libration  points  is  always  unstable  due  to 
the  existence  of  a  positive  real  root  of  the  characteristic  equation  for  any 
value  of  /i.  An  unperturbed  object  in  an  orbit  around  a  collinear  libration 
point  will  move  away  from  that  point. 

The  four  initial  conditions  corresponding  to  the  in-plane  motions  (X-Y 
plane)  cannot  be  arbitrarily  selected.  This  will  be  further  explained  in  Sec¬ 
tion  5.  Szebehely  discusses  that  the  selection  of  initial  conditions  at  the 
collinear  libration  points  show  considerable  inherent  instability.  Even  when 
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Figure  3:  Relative  Position  with  respect  to  L2. 

initial  conditions  are  chosen  according  to  the  linearized  requirements,  the 
spacecraft  orbits  maybe  once  or  twice  around  L2  before  moving  away.  Fortu¬ 
nately,  periodic  orbits  may  be  obtained  by  slight  modification  of  the  initial 
conditions  which  follow  from  the  linearized  solution.  The  point  is  that  the 
orbits  are  very  sensitive  to  changes  in  the  initial  conditions. 

4  Full  Nonlinear  Equations  for  Circular  Re¬ 
stricted  Three-Body  Problem  (relative  to 
a  libration  point) 

We  continue  from  Hamilton. 

Motion  around  a  collinear  libration  point  is  more  easily  expressed  in  a 
local  coordinate  frame,  shown  in  Figure  3,  with  the  origin  located  at  the 
libration  point  of  interest. 

Note  the  vector  addition 

f  —  f0  +  p 

where  fo  has  components  (A"o,  Yq,  Zq)  to  define  the  location  of  the  libration 
point  relative  to  the  system  barycenter.  The  “0”  subscript  refers  to  any 
one  of  the  libration  points.  (To  and  Z0  are  both  zero  for  collinear  libration 
points.)  Vector  p  has  components  (X,  Y,  Z)  to  specify  the  location  of  m3 
relative  to  the  libration  point.  Using 

X  =  X0  +  x  (4) 
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Y 

Z 


Y0  +  y 

Z0  +  z, 


(5) 

(6) 


be  reminded  that  (X,  Y,  Z)  is  the  spacecraft  position  relative  to  the  system 
barycenter.  The  radial  direction,  x,  is  collinear  with  ai,  going  from  the  Sun 
(or  largest  body)  through  the  libration  point.  The  cross-track  direction,  z, 
is  parallel  and  in  the  same  direction  as  a3,  the  orbit  normal.  The  component 
y  is  in  the  direction  of  a2 •  When  circular  motion  is  assumed,  the  y  direction 
is  the  tangential  velocity  direction  of  the  libration  point  around  the  system 
barycenter.  These  axes  form  an  orthogonal  coordinate  system.  Motion  in 
the  x-y  plane  is  referred  to  as  in-plane  and  any  motion  in  the  z  direction  will 
be  referred  to  as  out-of-plane. 

The  equations  below  describe  the  spacecraft’s  motion  in  a  local  coordinate 
frame  relative  to  the  libration  point  (origin)  chosen  by  selection  of  (X0,  Yq, 
Z0).  Unique  distances  corresponding  to  each  of  the  five  libration  points  exist 
for  each  set  of  (X0,  Y0,  Z0). 

The  equations  below  are  the  full,  coupled,  nonlinear  equations  of  motion 
for  the  circular  restricted  three-body  problem.  The  coordinate  frame  rotates 
with  a  constant  angular  speed. 

Equations  (1)  -  (3)  yield  by  simple  substitution  of  Equations  (4)  -  (6): 


x  —  2uy  —  cu2(X0  +  x)  = 
y  +  2ux  -  u2(Yv  +  y)  = 


H  i(Xq  +  X  +  D\)  p2  (Xq  +  X  —  -D2) 


r3 

'  1 


,r3 

'  2 


l^i(Yo  +  y)  y2(Yo  +  y) 


r.3 

'1 


r3 

'  2 


Z  = 


yi(Z0  +  z)  y2(Z0  +  z) 


where  rq  and  r2  are  now  defined  below  as 


f  1  =  \f  (X0  +  x  +  Drf  +  (To  +  y)2  +  (Z0  +  z)2 
r2  =  y/(X0  +  x~  D2)2  +  (Yq  +  y)2  +  (Z0  +  ^)2  . 
The  initial  conditions  to  be  specified  are 


x(0)  i(0) 

y( o)  2/(0) 
z(0)  i(0). 
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These  equations  are  derived  in  Appendix  B. 

For  the  L2  point,  where  only  Xq  ^  0,  the  equations  further  simplify. 


5  Linear  Equations  for  Circular  Restricted 
Three-Body  Problem 

At  a  chosen  equilibrium  point,  the  nonlinear,  full,  equations  of  motion  can  be 
linearized  in  order  to  exploit  standard  solutions  to  linear  equations.  The  lin¬ 
ear  equations  of  motion  about  the  local  coordinate  frame  (origin  at  collinear 
libration  point  L2)  are  given  here: 


x  —  2c vy  —  Uxxx  =  0 
y  +  2ux  —  UyyV  =  0 
z  -  Uzzz  =  0, 

where,  evaluated  at  the  L2  collinear  libration  point,  where  X  =  Xq,  Y  =  0, 
Z  =  0 

2  2^i  2/i2 

^  (X0  +  D ,)3  +  (X0  -  D2f 

2 _ hi _ h2 

1  {X0  +  D\)3  (. X0  -  D2 )3 

hi  h  2 

“(Xo  +  Z^i)3  {XQ-D2f 

The  initial  conditions  to  be  specified  are 

x(0)  x(0) 

</(o)  m 

2(0)  2(0). 

These  equations  are  derived  in  Appendix  C. 

When  evaluated  at  L-2  Uxx,  Uyy ,  and  Uzz,  which  were  formed  from  U, 
are  constant.  The  pseudopotential  U  is  the  centrifugal  plus  gravitational 
force  potential  defined  as 

U  =  1uj\X2  +  Y'2)  +  ^  ^ 

2  r  i  r2 


Uxx\L2  = 
uyy\l2  = 
Uzz\l2  = 
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and 


ri  =  yJ(X  +  D\)2  +  Y2  +  Z2 
r2  =  ^/(X  -  Z)2)2  +  Y2  +  Z2. 

The  respective  double- lettered  subscript  indicates  the  second  derivative 
of  U  with  respect  to  that  lettered-variable. 


The  selection  of  initial  conditions  for  the  first-order  linear  model 

summarizes  the  discussions  of  Szebehely  [3],  Farquhar  [5],  and  Wie  [4],  The 
in-plane  characteristic  equation  of  the  linearized  equations  of  motion  about 
the  collinear  equilibrium  points  is  of  fourth  degree.  Solution  yields  two  real 
roots  (equal  in  magnitude,  but  opposite  in  sign)  and  two  imaginary  roots. 
The  in-plane  position  has  a  convergent  mode  (due  to  the  negative  real  root), 
a  divergent  mode  (due  to  the  positive  real  root),  and  an  oscillatory  mode. 
All  three  collinear  libration  points  are  unstable. 

The  out-of-plane  characteristic  equation  of  the  linearized  equations  of 
motion  about  the  collinear  equilibrium  points  is  of  second  degree.  The  two 
eigenvalues  are  imaginary.  Thus,  out-of-plane  motion  is  simply  oscillatory. 

For  the  linearized  equations  of  motion,  these  authors  show  that  a  solution 
to  the  equations  can  be  made  to  contain  only  the  oscillatory  modes  with 
proper  choice  of  initial  conditions: 

Choose 


i(0)  =  \k)m 

(7) 

y(  o)  =  —ku!Xyx(0) 

(8) 

where  k  =  uix^U^x  and  uxy  is  the  nondimensional  frequency  of  the  in-plane 
oscillatory  mode  calculated  from 

^xy  =  y  Pi  +  \jp\  +  Pi 
Pi  =  2  —  (Uxx  +  Uyy)/2 
Pi  =  ~ UxxUyy , 
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which  means  that  once  initial  conditions  x(0)  and  y(0)  are  selected,  the 
corresponding  initial  velocity  components  cannot  be  chosen  at  will. 

Although  the  choice  of  the  z  initial  conditions  is  arbitrary,  Wie  states 
that  choosing  a:(0)  =  ;j(0)  =  0  and  i(0)  =  —y{Q)ujz  (with  cuz  =  \J\Uzz\),  the 
solution  further  reduces  to  a  quasi-periodic  Lissajous  trajectory.  We  will  not 
discuss  this  trajectory  at  this  time.  We  just  list  his  example  initial  condition. 

Szebehcly  continues  by  noting  Equations  (7)  -  (8)  do  not  hold  when 
higher-order  terms  in  the  function  U  are  retained.  Nevertheless,  the  set  of 
initial  conditions  given  by  these  equations  furnish  the  starting  point  of  an 
iteration  which  leads  to  the  proper  initial  condition  for  the  nonlinear  case. 
The  closer  the  initial  point  (x(0),  2/(0))  is  to  the  chosen  libration  point, 
the  better  a  differential  correction  scheme  will  work  to  furnish  the  initial 
conditions  desired  for  the  establishment  of  periodic  orbits  in  the  nonlinear 
problem. 

The  linearization  process  allows  only  first-order  terms  in  the  coordinates 
and  velocities;  consequently,  meaningful  results  must  be  connected  with  small 
values  of  these  variables.  Even  when  the  initial  conditions  are  chosen  accord¬ 
ing  to  the  linearized  requirements,  the  spacecraft  perform  only  one  or  two 
orbits  before  they  depart  from  the  area.  Again,  periodic  orbits  may  be  ob¬ 
tained  by  slight  modification  of  the  initial  conditions  which  follow  from  the 
linearized  solution.  The  point  is  that  the  orbits  are  very  sensitive  to  changes 
in  the  initial  conditions. 

The  linear  model  is  amenable  to  linear  analysis  methods,  which  are  more 
tractable  than  nonlinear  solutions.  The  linear  model  is  most  helpful  for 
preliminary  control  analysis  as  discussed  by  Hamilton.  The  nonlinear  insta¬ 
bilities  of  this  system  must  be  restrained  to  maintain  the  desired  periodic 
orbit  motion. 


6  Quadratic  Equations  for  Circular  Restricted 
Three-Body  Problem 

We  went  one-step  beyond  the  linear  dynamics  model  of  Section  5.  By  further 
expanding  U,  we  obtained  a  second-order  model  for  the  restricted  problem. 
This  model  should  yield  trajectories  closer  to  the  nonlinear  model  than  does 
the  linear  model. 
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where 


x  —  2  uij  —  u j2x 

y  +  2ux  —  u2y 

z 


2  Ax  —  ^B(  2x2  —  y2  —  z2) 

— Ay  +  3  Bxy 
— Az  +  3  Bxz 


A 

B 


y  i  ih 

(X0  +  D,) 3  +  (X0  -  D2f 

hi  h2 

(X0  +  L>i)4  +  (X0  -  -D2)4 


The  initial  conditions  to  be  specified  are 


x(0)  i(0) 

</(o)  m 

z(0)  i( 0). 


These  equations  are  derived  in  Appendix  D. 


7  Linear  Relative  Motion  Equations  for  Cir¬ 
cular  Restricted  Three-Body  Problem 

Consider  the  linearized  equations,  with  initial  conditions  selected  such  as  to 
avoid  exciting  the  divergent  mode.  The  solution  then  takes  the  form 


X 

x0  cos  fl(t  —  t0)  +  ^  sin  Q(t  —  1 0) 

r  = 

y 

= 

— kx0s\n£l(t  —  to)  +  y0cos£}(t  —  t0) 

z 

zq  cos u(t  —  t0)  +  —  sina;(t  —  to) 

where  subscript  0  refers  to  the  value  at  time  t  —  t0.  This  solution  can  be 
written  as  a  pair  of  decoupled  matrix  equations  of  similar  form: 


X 

cosf2(t  — t0)  |sinh2(t  — t0) 

X0 

.  y . 

— A:  sin  0(t  —  t0)  cosf2(t  — t0) 

.  y° . 

z 

cosa;(t  —  t0)  ^  sino;(t  —  to) 

zo 

z 

—uj  sin u> (t  —  t0)  coso;(t  — t0) 

.  . 
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For  notational  convenience,  let  the  vectors  x  and  z  refer,  respectively,  to 
[x  y]T  and  [z  i]T.  Therefore,  Equations  (9)  and  (10)  may  be  written  as 

x  =  -4fc,n(Mo)xo  (11) 

z  =  AU):UJ(t,to)z0,  (12) 

where 

T  cosfi(t-t0)  |  sin  Q(t  —  t0) 
k'n  ’  0  —  ksm£l(t  —  to)  cosfl(t  — to) 

is  the  generic  state  transition  matrix.  Because  of  the  obvious  similarity,  the 
relative  motion  solution  may  be  developed  for  one  of  these  matrix  equations, 
and  the  result  applied  to  the  other. 

7.1  Relative  Motion  of  Two  Objects  on  Same  Trajec¬ 
tory 

Now,  say  that  there  are  two  objects  which  satisfy  this  set  of  equations.  Object 
1  has  the  initial  conditions  which  have  been  discussed,  but  object  2  trails 
object  1  by  a  time  delay  At.  Therefore, 

r 2(t)  =  ri(t  -  At), 

with  the  subscript  referring  to  the  object  number.  Let  the  relative  motion 
be  given  by 

Ar  =  ri  -  r2, 

with  the  initial  state 

ri  (t0)  =  1*2  (t0  +  At)  =  r0. 

First  consider  the  in-plane  motion  of  Equation  (11).  The  in-plane  motion 
of  objects  1  and  2  is  given  by 

xi  =  Akjn(t,  t0)x0 
x2  =  Ak>n(t  —  At,  t0)x0. 

Therefore,  the  relative  motion  is 

Ax  =  [Akjn(t,  t0)  -  ^4fc,n(t  -  At,  t0)]x0.  (13) 
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Now, 


Ax0  =  Ax(t0)  =  xi0  -  x20 

=  Ak,n(to,  t0)x0  —  Ak^(t0  —  At,  fo)xo 
=  [I  -  Ak,n(t0  -  At,  t0)]x0. 

And  so, 

x0  =  [I  -  Akin(t0  -  At,  f0)]_1Ax0.  (14) 

Substituting  into  Equation  (13), 

Ax  =  [Ak^(t,  to)  —  Ak)n(t  —  At,  to)]  [I  —  Ak,n(to  —  At,  to)]  1Ax0. 

From  the  properties  of  the  state  transition  matrix, 

Ak,n(t  —  At,  t0)  =  Akfi(t,  t0)Ak^{t0  —  A t,t0). 

Therefore, 

Ax  =  [Ak^{t,  to)  —  Ak£i(t,  to)Ak^{to  —  At,  t0)] 

[I  ~  Akjn(to  —  At,  to)]  1Ax0 

=  Ak,n(t,  t0)[/  —  Akjci(to  —  At,  t0)][/  —  Ak,n(to  —  At,  to)]  xAx0 
=  -Afc,n(t,  t0)  Ax0. 

The  same  result  may  be  applied  to  Equation  (12).  The  resulting  relative 
motion  equations  are  then 

Ax  =  Afc;n(t,t0)  Ax0 
Az  =  A»^(t,to)  Az0, 


with  the  full  matrices  as  given  in  Equations  (9)  and  (10). 

It  may  be  preferable  to  present  the  relative  motion  in  terms  of  the  relative 
position  at  the  time  of  insertion  of  the  second  object  —  that  is,  at  time  t0+At. 
Referring  to  Equation  (13), 

Ax(t0  +  At)  =  [Ak)n(t0  +  At,  t0)  -  Ak^(t0,  t0)]x0 
—  [Ak,n(to  +  At,  t0)  —  7]x0. 
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Following  the  same  logic  as  above, 

Ax0  =  [Afc)n(to  +  At,  t0)  -  /]_1Ax(f0  +  At), 

and  so 

Ax  =  [Ak^(t,  t0)  -  Akfl(t  -  At,  t0)]x0 

=  [-4fc,n(t  —  At,  to)^4fc,n(to  +  At,  to)  —  Akjn(t  —  At,  to)] 

[Ak,n(to  +  At,  to)  —  I]  1Ax(t0  +  At) 

=  Akjn(t  —  At,  t0)Ax(t0  +  At). 

Once  again,  the  same  result  applies  to  the  out-of-plane  motion,  giving 

Az  =  A^(t  -  At,  t0)Az(t0  +  At). 

Note  that  these  results  could  also  be  written  as 

Ax  =  Ak>n(t,  t0  +  At)  Ax(t0  +  At) 

Az  =  A^it,  t0  +  At)Az(t0  +  At). 

7.2  Relative  Motion  With  Insertion  Position  or  Time 
Error 

Next,  consider  the  possibility  of  an  injection  error  for  the  second  object,  such 
that  it  does  not  in  fact  have  the  same  initial  state  as  the  first.  Let  X05  refer 
to  the  in-plane  state  error  at  injection,  such  that 

x2(f0  +  At)  =x0  +  x06. 

The  resulting  relative  state  is  then,  simply, 

Ax  =  Akfl(t,t0  +  At)(Ax(t0  +  A t)a  -  xoe,), 

where  A x(t0  +  A t)a  is  the  nominal  relative  state  at  the  second  object’s  in¬ 
jection  time,  as  discussed  above.  Again,  the  same  result  applies  to  the  out- 
of-plane  motion,  giving 

Az  =  AU)U)(t,to  +  At)(Az(to  +  A  t)a  +  z0b). 
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It  may  also  be  useful  to  examine  an  error  in  the  injection  time  of  the 
second  object.  Let  the  time  error  be  A tb,  such  that 

At  =  A  ta  +  A  tb, 

where  A ta  refers  to  the  nominal  delay  between  the  insertion  of  the  first  and 
second  objects.  In  this  case,  the  in-plane  state  of  the  second  object  is 

X2  =  Akp(t  -  At a  -  A tb). 


The  relative  state  is  then 


Ax  =  [Akp(t,  t0)  -  Akin(t  -  A ta  -  A tb,  t0)]x0.  (15) 

At  the  nominal  injection  time,  to  +  A ta,  the  relative  state  is 

Ax(t0  +  A ta)  =  [Ak,n(t0  +  Ata,  t0)  -  Ak^(t0  -  A tb,  t0)]x0.  (16) 

The  second  state  transition  matrix  may  be  expanded  in  a  Maclaurin  series 
about  A  tb  =  0.  A  first-order  expansion  gives 


Ak^(t0  —  A  tb,  to) 


1  + 


0  —  £A  tb 

kftAtb  0 


=  I  +  Bk,n(Atb). 


Substituting  into  Equation  (16), 


Ax(t0  +  Ata)  =  [Afc,o(t0  +  Ata,  to)  -  L]x0  -  5fcjn(Atfe)x0. 


Note  that  the  first  term  gives  the  nominal  relative  position  at  the  nominal 
insertion  time  of  the  second  object;  the  second  term  gives  the  error.  Denote 
this  first  term  as  A x(f0  +  A ta)a-  Then, 


x0  =  [Ak,si(to  +  Ata,  t0)  -  I]  1A^(t0  +  A ta)a. 

Substituting  into  Equation  (15),  and  again  invoking  the  properties  of  the 
state  transition  matrix, 

Ax  =  [Aktn(t,  t0)  -  AkJ[i(t  -  A ta  -  A tb,  t0)] 

[Aktn(to  +  Ata,  to)  —  I]  1Ax(t0  +  A ta)a 
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—  [A k,n(t,  to)  —  Ak,n(t  —  A ta,  to)^4fc,n(to  —  At&,  to)] 

[^fc,n(to  +  A  ta,  to)  -  /]_1Ax(t0  +  Ata)a 
=  {Ak,n(t,  t0)  —  Ak^(t  —  A ta,  t0)[/  +  -Bfc,n(Atft)]} 

[^fc,n(to  +  A ta,  to)  -  /]_1Ax(t0  +  A ta)a 
=  Ak:n(t  —  A ta,  to){Afc;n(to  +  Ata,  t0)  —  [/  +  -Bfe,n(At)]} 

[Afc,n(to  +  Ata,  to)  —  /]  1Ax(t0  +  A  ta)a 
=  -4fc,n(t  —  Ata,  t0){/  —  -Bfc)n(At)[AfcjQ(t0  +  Ata,  t0)  —  /]  1} 
Ax(t0  +  Ata)a 

=  Afc;n(t  -  A ta,  t0)  Ax(t0  +  Ata)a 

—  Ak,n(t  —  Ata,  t0)-Bfc,n(At)[Afcin(t0  +  A  ta,t0)  —  /]  1 
Ax(t0  +  Ata)  a- 

Here,  note  that  the  first  term  is  again  the  nominal  relative  position,  while 
the  second  term  is  the  resulting  error. 

8  Some  Examples 

Currently,  there  are  four  different  MATLAB  scripts  that  present  numerical 
solutions  to  these  models: 

•  full  equations  of  motion  about  L2 

•  linear  equations  of  motion  about  L2 

•  quadratic  equations  of  motion  about  L2 

•  relative  motion  of  two  objects  on  the  same  trajectory. 

We  first  duplicated  Hamilton’s  results  shown  in  his  Section  8.1.1  to  check 
the  linear  code. 

For  subsequent  analysis,  we  used  a  consistent  set  physical  quantities  taken 
from  Dunham  and  Muhonen  [6]. 

Let  us  examine  the  same  initial  conditions  applied  to  these  four  numerical 
solutions.  We  begin  near  the  L2  equilibrium  point.  Consider  displacing 
the  spacecraft  10,000  km  away  from  L2  along  only  the  x  direction,  with 
7/(0)  =  —1,127.8  krn/day  as  defined  from  Equation  (8)  of  Section  5.  All 
of  y,  i(0),  z,  and  i(0)  are  zero.  As  noted  earlier,  motions  in  the  Z-plane 
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x  position  (km) 

- Full  - Linear - Quadratic 


Figure  4:  Solutions  to  Models  Show  Contrasting  Fidelity. 


are  uncoupled  from  the  X-Y  plane;  therefore,  for  these  initial  conditions,  all 
motions  will  be  in  the  X-Y  plane. 

Figure  4  contrasts  the  differences  among  the  three  models  of  the  equations 
of  motion  as  shown  for  170  days  (nearly  a  full  period  of  178  days).  The  full 
model  does  not  even  close  on  itself.  For  these  intitial  conditions  the  full 
model  diverges  and  the  quadratic  model  shows  very  good  agreement  with 
the  full  model.  The  linear  model  does  close  on  itself,  but  only  represents  the 
actual  motion  for  about  one-quarter  of  an  orbit.  Shown  on  the  figure  is  the 
L2  point,  which  is  the  origin,  and  the  starting  location.  MATLAB’s  “ode45” 
solver  was  used,  which  is  based  on  an  explicit  Runge-Kutta  (4,5)  formula. 
The  print  interval  was  one  day. 

Continuing  with  the  same  data  produced  for  Figure  4,  we  turn  to  Figure  5. 
Here  we  take  the  component  positions  of  the  linear  and  quadratic  models  and 
subtract  the  corresponding  position  of  the  full  model.  Although  plotted  for 
22  days  beyond  one  period,  we  can  see  how  rapidly  the  divergence  grew  as 
the  second  orbit  began. 
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Figure  5:  Relative  Errors  Among  Models  Show  Significant  Position  Differ¬ 
ences. 
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Figure  6:  Relative  Position  of  Two  Spacecraft. 


We  move  on  to  the  analytic  solution  of  the  relative  motion  between  two 
spacecraft.  Here  we  used  the  linear  solution  coded  from  Equation  (13)  of 
Section  7.1.  See  Figure  6,  which  shows  the  motion  of  spacecraft  ^2  with 
respect  to  spacecraft  #1.  The  calculations  based  on  the  initial  conditions 
discussed  above  were  made  for  each  spacecraft.  They  have  the  same  initial 
position  and  velocity  conditions  yet  spacecraft  ^2  began  one  day  after  space¬ 
craft  $T.  Each  diamond  on  the  figure  shows  the  position  of  spacecraft  ^2 
in  steps  of  two  days.  You  might  consider  yourself  to  be  at  the  origin  riding 
on  spacecraft  #1.  With  a  one-day  delay,  the  position  of  spacecraft  ^2  is 
plotted  every  two  days  (to  offer  clarity  between  plot  symbols).  At  the  orbit’s 
beginning  and  halfway  points,  the  relative  change  of  position  of  spacecraft 
7^2  is  less  than  other  orbit  points. 
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Appendices 


A  Derivation  of  the  Full  NonLinear  Equa¬ 
tions  of  Motion  for  the  Circular  Restricted 
Three-Body  Problem  (relative  to  barycen- 
ter) 

The  development  of  the  equations  of  motion  will  follow  that  of  Hamilton.  In 
our  reading  of  his  thesis  to  understand  his  work,  we  rederived  the  equations 
of  motion.  More  details  are  given  in  these  respective  appendixes  to  maintain 
a  record  of  our  labors. 

Returning  to  Figure  1,  which  is  a  rotating  coordinate  system,  we  find  the 
need  to  add  an  inertial  reference  frame  to  define  the  rotation  vector. 

An  inertial  coordinate  system  is  defined  with  origin  at  the  system  barycen- 
ter.  The  orthogonal  basis  vectors  n \ ,  77.2 ,  and  77,3  are  fixed.  Vectors  h\  and 

77.2  lie  in  the  same  plane  as  the  rotating  basis  vectors  aq  and  a2-  The  vector 

77.3  is  always  aligned  with  03.  Once  every  revolution,  fq  is  aligned  with  cq 
at  the  same  time  that  77,2  is  aligned  with  a2.  With  respect  to  the  inertial 
coordinate  frame,  the  two  large  bodies  are  rotating  about  their  barycenter 
with  a  constant  angular  velocity 


U) 


0 

0 


n 


U) 


cc m3 


where 

G(m1  +  m2) 

U  ~  V  D 3  ’ 

and  G  is  the  gravitational  constant. 

With  rotating  reference  axes,  the  calculation  of  the  spacecraft’s  velocity 
with  respect  to  an  inertial  frame  is 

dr  dr 

—  —  —  +  00  x  r 

Clt  n  Cit  a 
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where  the  position  vector  of  the  spacecraft  in  the  rotating  frame  is 


T  —  A  d\  A  Y  S  2  A  Zd  3. 


Substitute  into  the  calculation  and  collect  terms: 


dr 

dt  n 


A  Si  F  Y  CL 2  +  ZS3  + 


0-2  ^2 
0  0  a; 

A  y  z 


r  —  (A  —  ycu)Si  A  (i  A  Acu)S2  A  ZS3. 


With  rotating  reference  axes,  the  calculation  of  the  acceleration  of  vector 
f  with  respect  to  an  inertial  frame  is 


d2r  ^  ,  _  .  dr  d2r 

—  =  u  x  rn  A  w  x  (u  x  ra)  A  2w  x  —  +  — 

(XL  n  (XL  a  (XL  a 


The  first  term  is  zero  because  the  system  rotates  at  constant  angular 
speed. 

Substitution  into  the  equation  yields 


d2r 

dt 2  n 


co  x  ( — c oY Si  A  00X0,2)  A  2 to  x  (A Si  A  h  cl 2  A  ZS3) 


A  (ASi  A  Y  S2  A  ZS3) 


CL\ 

S2 

S3 

Si 

S2 

S3 

0 

0 

00 

A 

0 

0 

2c jo 

A  (A  Si  A  Y S2  A  ZS3 

-uY 

coX 

0 

A 

Y 

Z 

f  =  (A  -  2cW  -  cu2A)Si  A  (y  A  2u;A  -  cu2y)S2  A  Za3 


We  associate  the  kinematics  of  f  with  the  force  of  gravity  acting  on  the 
system.  Gravity  is  the  only  force  acting  on  the  system. 

From  Newton’s  second  law,  which  states  that  the  time  rate  of  change  of 
linear  momentum  is  proportional  to  the  force  applied.  For  our  fixed-mass 
system  of  interest 


F  =  m:ir  = 


Grn  1  rrcs  _ 
- r 

N3  - 


Grri2m3 


T2- 
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As  a  reminder,  rn  \  is  the  mass  of  the  larger  large  body,  m2  is  the  mass  of 
the  smaller  large  body,  and  m3  is  the  mass  of  the  spacecraft.  Now  substitute 
into  this  expression  the  equations  for  f,  fq  and  r2  and  separate  the  vector 
equation  into  the  three  components. 


(X  -  2 ujY  -  u2X)di  +  (Y  +  2ujX  -  uj2Y)a2  +  Za3  = 

— ; — or  [(A  A  Di)di  +  Ya2  +  Za, 3]  — , — j^r2[(X  —  D2)a  1  +  Y cl2  +  Zci3] 
\t  ir  p  2 1 3 


X-2  uY-u2X  = 
Y+2uX-u2Y  = 
Z  = 


fii(X  +  D\)  n  2  (A  —  D2) 


|n|3 

A^i  Y  n2Y 


r  2 


lrl|: 


1"2 


/i  1  Z  fi2Z 


r  2 


From  Wie  [4]  we  note  these  equations  of  motion  can  also  be  expressed  in 
terms  of  a  pseudopotential  U  =  U(X,  Y,  Z)  as  follows: 

X-2  uY  =  %  =  l 'Jx 
oX 

Y  +  2uX  =  %  =  UY 
01 

7  =  dU  a  n 

dz  Z: 

where  the  pseudopotential  U,  which  is  the  centrifugal  plus  gravitational  force 
potential,  is  defined  as 

u  =  in*2  +  W)  +  -  +  -■ 


n  r2 


Also 


n  =  \J  (X  +  Di)2  +  Y2  +  z2 


r2  =  y/(X-D2)*  +  Y2  +  Z*. 
Again,  the  initial  conditions  to  be  specified  are 


WO)  WO) 
WO)  Wo) 
WO)  Wo). 


26 


B  Derivation  of  the  Full  Nonlinear  Equations 
of  Motion  for  the  Circular  Restricted  Three- 
Body  Problem  Near  L2 

Begin  with  the  results  of  the  derivation  of  the  full  nonlinear  equations  of 
motion  relative  to  barycenter,  Equations  (1) — (3): 


X-2  ujY-uj2X  = 
Y  +  2uX-uj2Y  = 
Z  = 

Substitute  using  (4)-(6): 

X 

Y 

Z 


m(x  +  p  i)  h2{x-d2) 

ki|3  |r2|3 

\mY_  _  IxY 

\r i|3  | r2 | 3 
Yi-Z  1^2  Z 

In  I3  V2 13 

=  X0  +  x 
=  Y0  +  y 
—  Zq  +  z, 


to  obtain 

(A"0  +  x)  —  2u  (Yq  +  y)  —  uj2(X0  +  x)  = 
(h)  +  y)  —  2uj(Xq  +  x)  —  o;2(yo  +  y)  — 

(Zq  +  z)  = 


l*i(x0  +  x  +  Di)  /i2(A  0  +  x  —  D2) 


! r  1 1 3  N3 

fii(Y0  +  y)  y2{Y0  +  y) 

13 


r  2 


Ah(-^o  +  z)  y2(Z0Zz) 


ln|J 


I  ^2 1 : 


By  definition  the  equilibrium  point  is  stationary;  therefore,  many  terms 
are  zero. 

A0  —  Ox  —  0 

Y0  =  0  y  =  0 

Z0  =  0  z  =  0. 

The  equations  reduce  to  the  following: 

//i(*o  +  x  +  Di)  ^(Xq  +  x  —  D2) 


x  —  2uy  —  ar(A0  +  x)  = 


,r3 

'1 


r3 

'  2 
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,  O  •  2  (\r  i  \  MYo  +  y)  ^(Y0  +  y) 

y +  2ux -u  (Y0  +  y)  = - 3 - 3 - 

rl  r2 

„  _  fli(Z0  +  z)  ^{Zq  +  z) 

Z  ~  Ti  Ti  ’ 

'  1  '2 

where  r\  and  r-2  are  now  defined  below  as 

r\  =  (A"0  +  x  +  -Di)2  +  (Vo  +  y)2  +  (Zo  +  z )2 

r2  =  7(X0  +  x  -  D2)2  +  (ho  +  y)2  +  (Zo  +  z)2. 

The  initial  conditions  to  be  specified  are 

x(0)  i(0) 

</(o)  m 
2(0)  2(0). 


C  Derivation  of  the  Linear  Equations  of  Mo¬ 
tion  for  the  Circular  Restricted  Three-Body 
Problem  Near  L2 

We  develop  linear  equations  to  have  models  applicable  for  linear  analysis 
methods. 

We  document  details  not  shown  in  Hamilton’s  thesis. 

We  repeat  the  nonlinear  equations  for  convenience. 


where 


Also 


X  -  2c uY 
Y  +  2  ojX 
Z 


dU 

dX 

dU 

dY 

dU 

~dZ 


U 


Y2(X2  +  Y 2)  +  ^ 
2  r  1 


+ 


l-h 

r2  ' 


ri  =  y/{X  +  D1)2  +  V2  +  Z2 
r2  =  y/{X  -  D2 )2  +  V2  +  Z2. 


We  substitute  the  individual  coordinates 


X  =  X0  +  x 
Y  =  Y0  +  y 
Z  =  Zq  +  z 


to  obtain 


(Wo  +  x)  —  2u(Yo  +  y) 
(Wo  +  y)  —  2u)(Xq  +  x ) 
(Zq  +  z) 


dU 

dX 

dU 

dY 

dU 

~dZ 


Ux 


UY 


Uz 
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Now,  perform  a  Taylor  series  expansion  on  these  partials  of  U,  about  the 
libration,  equilibrium,  point.  The  first  two  terms  of  this  series  for  a  function 
of  one  variable  are 

f(x)  =  f(X0)+xf'(X0). 

For  Ux,  let  f(x,y,z )  = 

Note  that  because  the  function  is  already  a  first  derivative,  we  stop  the 
Taylor  series  at  its  first  derivative  so  we  can  stop  the  series  at  the  quadratic 
term.  We  have  a  function  of  three  variables  so  the  Taylor  series  is  modified 
as  follows  (with  L2  substituting  for  subscript  0  as  the  chosen  expansion  and 
evaluation  point) 

d  d 

f(x,y,z)  =  Ux\L2  + x—Ux  +  y-^Ux 

L2 

For  UY ,  let  f(x,y,z )  =  (p,  then 

d  d 

f(x,  Vi  z)  =  Uy\L2  +  x-q^uy  +  ygYUY 

L2 

For  Uz,  let  f(x,y,z )  =  ||,  then 

d  d 

f(x,  V’ z)  =  uz\l2  +  xoxUz  +  y9YUz 

L2 

By  dehnition  at  the  equilibrium  point,  there  is  no  motion.  Many  terms 
are  zero. 

XL2  =  xL2  =  U x\l2  —  0 
Yl2  =  Vl2  =  Uy\l2  =  0 
Zl2  —  zl2  —  Uz\l2  —  0 

Additionally,  terms  of  Ux  and  Uy  involving  the  partial  of  Z  are  zero  be¬ 
cause  motions  in  the  Z-direction  are  uncoupled  from  the  X-Y  plane  motions. 

Therefore 

x  —  2uy  =  xUxx\l2  +  uUyx\l2 
V  -  2cni;  =  xUXy\l2  +  yUYY \l2 
z  =  zUzz\l2 


+  zlzUz 


d 

+  zazUr 


d 

+  ZdZUx 
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The  above  equations  are  the  first-order,  or  linear,  coupled,  equations  of 
motion  for  the  circular  restricted  three-body  problem.  These  describe  the 
motion  of  the  spacecraft  relative  to  an  equilibrium  point  (X0,Y0,  Z0).  The 
equilibrium  point  rotates  with  a  constant  angular  velocity. 

U\o  =  y(X02  +  y02) 


For  the  restricted  three-body  problem  at  L2,  the  Yq  and  Zq  terms  are  zero 
and  the  equations  reduce  to 

x-2uy  =  xUXx\l2 
y  +  2ux  =  ijUyy\l2 
z  =  zUzz\l2 
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Return  to  the  general  form  of  U  as  we  first  carry  out  the  derivatives  of 
U,  then  substitute  for  the  libration  point. 


U  =  y  (X2  +  X 2)  + 
For  X: 

dll 

dX 


hi 


+ 


h  2 


yJ(X  +  DJ2  +  Y2  +  Z2\  \y/(X  -  D2)2  +  Y2  +  Z2 


yJ(X-D2)*  +  Y2  +  Z\ 


=  u2X 


hi 


h  2 


0*7 
<9X  L 
at/ 
dX  L 


((X  +  A)2  +  x2  +  z2)-* 
({X  -  D2f  +  X2  +  Z2)-1* 


Employ  the  differential  d(un)  =  nun  1du  to  calculate  the  last  two  terms 
-he  equation  : 

ftH[((A'  +  £>,)2  +  l"2  +  Z2)-i 
=  +  Bj)2  +  Y-  +  Z2)-iX((X  +  D{f  +  Y2  +  Z2) 


(2X  +  2  Dx 


_hi  _ 1 _ 

2  yJ(X  +  Dj2  +  Y2  +  Z 2)  3 

fii(X  +  Di) 

^/(x  +  d i)2  +  x2  +  z2)  3 

Of  course  this  calculation  is  similar  for  the  term  beginning  with  [i 2.  Here 
assign 


n  =  ^(x  +  z^)2  +  x2  +  z2 

r2  =  J(X  -  X>2)2  +  X2  +  Z2 
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Finally 


dU=  2  MX-KPi)  fj,2(X-P2) 

dX  U  rf  r2 

Now  take  the  second  derivative  using  d(uv)  =  udv  +  vdu: 


=  a;2-/u1|(X  +  JD1)  [-^((X  +  JD1)2  +  F2  +  Z2)-f(2X  +  2JD1) 

+  ((X  +  D1)2  +  F2  +  Z2)-i(l)| 

-  /i2  |(X  -  F>2)  -  D2)2  +  F2  +  F2)-i(2X  -  2D2) 

+  ((X  -  D2)2  +  F2  +  Z2)-i(l)| . 


Now  reduce  the  terms  and  substitute  for  ry  and  r2 


U  XX  —  CO 


1  3{X  +  £>i); 


1  3(X  -  D2 


For  F: 

and 

For  F: 

and 


5f/=  2  hiF  /i2F 

<9F  U  rf  r\ 


dU  Hi  Z  ii2Z 
dZ  r?  r2 
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If  desired,  evaluate  these  at  a  collinear  libration  point,  say  L2  where 

X  =  xQ,  y  =  0,  z  =  0. 


Uxx\L2 

uyy\l2 

Uzz\l2 


L U2  + 


2hi 

(Xo  +  D 1)3 


00 


hi 

(Xo  +  Di)3 

hi 


(X0  +  A)3 


2h2 

'  +  (Xo  -  D2f 
h  2 

~  (X0  -  ZZ2)3 
h2 

(X0  -  D2)3 
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D  Derivation  of  Quadratic  Differential  Equa¬ 
tions  for  Circular  Restricted  Motion  Near 

L2 


We  repeat  the  nonlinear  equations  for  convenience: 


X  -  2c uY 
Y  +  2  ojX 
Z 


dU 

dX 

dU 

dY 

dU 

~dZ 


where 


U  =  %-{X2  +  Y2)  +V±  +  ^L 
2  r  i  r2 


and 


n2  =  (X  +  Dif  +  Y2  +  Z2 
r22  =  (X  -  D2)2  +  Y2  +  Z2 


(17) 

(18) 

(19) 

(20) 

(21) 

(22) 


The  right  sides  of  the  equations  of  motion  can  be  formed  from  differenti¬ 
ation  of  Equation  (20): 


dU  _  2  Ah  &n  n2  dr 2 

dX  U  ri2  dX  r22  dX 
dU  _  2  hi  dr i  ji2  dr2 

dY  ri2  dY  r22  dY 

dU  Hi  dr i  Hz  dr2 

dZ  =  ~V^~dZ  ~  V^dZ' 


(23) 

(24) 

(25) 


The  required  derivatives  of  r i  and  r2  are  formed  by  differentiating  Equa¬ 
tions  (21)  and  (22): 

dri  _  X  +  D1 

dX  r'i 

dr  i  Y 

dY  r\ 

dr  i  Z 

dZ  ri 
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dr2  =  X  -  P2 
dX  r2 

dr  2  Y 

d  Y  r2 

dr  2  Z 

dZ  r2 

Equations  (23)  -  (25)  are  expanded  in  Taylor  series  expansions  about  the 
L2  libration  point,  at  coordinates  (X,Y,Z)  =  (X0,0,  0).  Let 


X 

1 

1 

0* 

_ 1 

y 

= 

Y 

z 

Z 

(26) 


give  the  displacement  from  L2.  Denote  dU/dX  by  Ux,  and  similarly  for  Uy 
and  Uz-  Then,  the  Taylor  series  through  second-order  in  x,  y,  z  are  given  by 


Ux  =  UX\L  + 


(d  Ux 

,  dUx 

d  Ux 

[  dx 

L  ay 

U/2 

y  + 

Li 

dZ 

1 

2 


d2U 


x 


dx 2 

+ 


X 2  + 


d2U 


x 


l2 

d2Ux 


dXdY 


dY 2 
xy  + 


y2  + 


d2U 


x 


l2 

d2Ux 


l2 


dXdZ 


dZ 2 
xz  + 


l2  / 

z2 

l2 

d2Ux 


l2 


Uy  =  Uy\l  T 


(dUy 

dUy 

dUy 

\  dx 

Lx+  BY 

U'2 

L  V+  dZ 
^2 

1 

2 


d2U 


5" 


dx2 
+ 


x2  + 


d2U 


Y 


L2 

d2Uy 


dY2 


y2  + 


d2U * 


y 


dXdY 


xy  + 


L2 

d2Uy 


l2 


dXdZ 


dZ2 
xz  + 


dYdZ 


Li  / 

L2 

d2Uy 


yz 


l2 


L2 


Uz  =  Uz. 


1  (duz 

dUz 

dUz 

\Li  +  ^  dx 

1 

2 


(d2Uz 

2  1 

d2Uz 

9  . 

d2Uz 

\  dX2 

x  + 
l2 

dY2 

to 

t 

+ 

dZ2 

+ 


d2U? 


dXdY 


xy  + 


d2U7 


l2 


dXdZ 


xz  + 


dYdZ 


l2  / 

^2 

L2 

d2Uz 


yz 


l2 


L2 


dYdZ 


yz 


l2 


(27) 


(28) 


(29) 
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A  vertical  bar  with  subscript  L2  indicates  evaluation  at  the  coordinates  of 
the  libration  point. 

The  derivatives  of  U  are  constructed  by  differentiating  Equations  (23)  - 
(25).  Note  that  higher-order  derivatives  of  rq  and  r2  are  required.  These  are 
given  by 


d  2rq 

1 

(x  +  do2 

dX 2 

n 

rq3 

<92rq 

1 

Y2 

dY2 

r'i 

r  i3 

d2r\ 

1 

Z2 

d  Z2 

rq3 

d  2rq  _  (X  +  DJY 

dXd  Y  ~  rq3 

d2rl  (. X  +  Di)Z 

dXdZ  ~  7? 

<92rq  YZ 


dZdZ  n3 

d 2r2  1  (X  -  D2 )2 


dX2 

r2 

r23 

d  2r2 

1 

Y2 

dY2 

T  2 

r23 

d2r2 

1 

Z2 

dZ 2 

r2 

r23 

d2r2 

(X 

-  t>2)e 

OXdY  r23 

d2r2  {X  -  D2)Z 


dXdZ 

d2r2  YZ 


dZdZ 


f  2 


3  ' 


The  evaluated  partial  derviatives  are  then: 


dUx 

dX 

dUY 

dY 


Li 

Li 


2  _| _ ^/i 1 _ I _ 2/i2 

^  +  (X0  +  £>i)3  +  (X0  -  D2)3 

2 _ hi _ h2 

'  "  (A'0  +  -Di)3  (A'o  -  Z)2)3 
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dUz 

hi 

h  2 

dZ  l2 

(X0+Dy 

)3 

* 

O 

1 

b 

to 

CO 

d2Ux 

6/I1 

6  h2 

dX2  T 

L  2 

(Xo  +  D, 

V4 

(X0  -  D2)4 

d2Ux 

3/ii 

I 

3/i'2 

dY2  T 
^2 

(X0  +  D 04 

1 

(X0  -  D2 )4 

d2Ux 

3/i  1 

1 

3/i  2 

dZ2  l2 

(X0  +  Di)4 

1 

(X0  -  D2 )4 

and 

dUx  _  dUx  _  d UY  _  d2Ux  _  d2Ux 
dY  T  ~  dZ  r  ~  dZ  T  ~  3XY  T  ~  dXZ  T 

L*2  L,  2  L  2  L  2  L  2 

d2Ux  _  d2UY  _  d2Uy  _  d2UY  _  d2Uz 

“  dYZ  ,  -  dY2  T  ~  dY Z  T  ~  dZ 2  T  ~  dZ 3  r  " 

L‘2  J^‘2  ^2  2 

Note  that,  because  of  the  continuity  of  the  derivatives,  any  other  derivatives 
are,  in  fact  repeats  of  those  listed.  For  example, 

d2UY  _  d2  (dU\ 
dX2  ~  dX^yW) 
d3U 
dX2Y 
d2  fdU\ 
dXY  V<9X ) 
d2Ux 
dXY ' 

By  definition,  the  L2  coordinates  satisfy  the  equilibrium  condition  of 
Equations  (17)  -  (19).  Therefore, 

Ux\l2  =  Uy \La  =  Ux\l2  =  0. 

Substituting  the  derivatives  into  Equations  (27)  -  (29), 


38 


UY 


Uz 


Define 


A4 1 


A4  2 


u; 


+  3 


[(X0  +  D1)4  + 

(X0  -D2)\ 

A4! 

A42 

(Xo  +  DO3 

(X0  -  D2f 

A4! 

A42 

(^0  +  A)4  (X0-D2)3 

A4!  A4 2 


(2a;2  -  y2  -  z2) 

y 

xy 


+  3 


(X0  +  Di)3  (X0  -  D2f 

A4!  .  A42 


+ 


(Xo  +  Di)4  (X0  -  D2f 


xz. 


A  k 
B  = 


A4 1  A42 

(X0  +  DO3  (X0  -  D2f 

/ii  A42 

(X0  +  D\)4  +  (X0  —  D2)4  ’ 


Then,  substituting  in  Equations  (30)  -  (32), 

Ux  =  (cn2  +  2A)x  —  2x2  —  y2  —  z2) 

Uy  =  (cu2  —  A)y  +  3  Bxy 
Uz  =  —  Az  +  3  Bxz. 


(30) 

(31) 


(32) 


(33) 

(34) 


Substituting  into  Equations  (17)  -  (18),  and  using  Equation  (26)  on  the  left 
side,  gives 

3 

x  —  2uy  —  uj2x  =  2Ax  —  -B{2x2  —  y2  —  z2) 

y  +  2ux  —  u2y  =  —Ay  +  3  Bxy 
z  =  —  Az  +  3  Bxz. 
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1  Introduction 


Some  continuing  long-term  goals  of  our  sponsor  are  to  develop  high-fidelity 
equations  of  motion  representing  the  formation  flying  of  spacecraft  near  the 
Sun-Earth  L2  point  and  equations  describing  the  relative  motion  between 
these  constellation  members.  The  equations  are  to  be  used  to  develop  orbit 
control  schemes. 

This  is  our  second  report,  which  continues  from  Part  1  [1],  with  further 
derivations  of  analytical  expressions  for  the  relative  motion  between  a  for¬ 
mation  of  spacecraft  orbiting  near  the  Sun-Earth  L2  point.  This  research  is 
loosely  motivated  by  formation  flying  concepts  for  the  MicroArcsecond  X-ray 
Imaging  Mission  (MAXIM)  Pathfinder. 

To  begin  the  understanding  of  the  behavior  of  objects  in  formation  near 
L2,  we  begin  with  the  assumptions  of  a  circular  restricted  three-body  prob¬ 
lem.  In  this  report,  we  are  not  modeling  the  true  eccentric  orbit  of  the  Earth 
and  Moon  about  their  respective  primaries  and  other  planetary  and  solar 
gravitational  orbit  perturbations. 

In  this  report,  we  present  a  preliminary  understanding  of  the  relative 
motion  between  two  spacecraft  orbiting  the  chosen  L2  point.  While  com¬ 
ponents  are  further  identified  below,  we  describe  the  motion  of  a  typical 
telescope  spacecraft  with  respect  to  the  central  hub  spacecraft  of  the  constel¬ 
lation.  The  hub  is  treated  as  the  constellation’s  reference  point;  its  orbital 
path  about  L-2  was  examined  in  Part  1.  The  description  of  motion  of  one 
telescope  spacecraft  with  respect  to  the  hub  spacecraft  can  be  applied  to  as 
many  telescope  spacecraft  as  needed. 

After  a  significant  literature  search,  we  state  that  what  we  are  uniquely 
contributing  is  a  description  of  spacecraft  in  formation  flight  about  the  Sun- 
Earth  L2  point. 

First,  we  list  the  components  of  the  MAXIM  mission.  Second,  we  develop 
the  differential  equations  of  motion  for  the  telescope  spacecraft  with  respect 
to  the  hub.  The  resulting  solution  is  developed  using  the  Lindstedt-Poincare 
perturbation  method  to  ensure  periodic  motion.  Periodic  motion  of  the  in¬ 
plane  and  out-of-plane  orbit  is  desirable  to  maintain  the  formation.  Finally, 
we  provide  analytic  solutions  of  relative  motion.  This  report  thus  contains 
the  following: 

•  development  of  differential  equations  of  motion  for  the  telescope  with 
respect  to  the  hub 
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•  description  of  the  linear  telescope  motion  about  the  hub 

•  discussion  of  the  preference  of  a  halo  orbit  (rather  than  a  Lissajous 
orbit)  of  the  hub  about  the  L2  point 

•  inclusion  of  linear  hub  motion  effects  upon  the  motion  of  the  telescope 
relative  to  the  hub 

•  inclusion  of  quadratic  hub  motion  effects  upon  the  relationships  be¬ 
tween  the  telescope  frequency  and  amplitude 


2  Problem  Definition 

The  current  configuration  of  the  MAXIM  Pathfinder  mission  is  depicted  and 
described  below. 


Figure  1:  Aperture  Formation  of  MAXIM  Pathfinder,  Concept  of  June  2002. 
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Seven  spacecraft  shown  in  Figure  1  are  in  the  same  plane  forming  a  flat 
aperture.  The  optical  hub  is  in  the  “center”  of  six  free-flying  spacecraft 
(telescopes).  Dashed  lines  do  not  indicate  a  physical  connection,  but  rather 
indicate  random  radial  directions  listed  ry  through  r§.  The  scalar  length  of 
r  ranges  from  100  to  500  meters.  Additionally,  the  six  free-flying  spacecraft 
are  loosely  spaced  60°  from  one  another. 

There  is  also  a  detector  spacecraft  located  approximately  20,000  km  and 
90°  out  of  the  plane  formed  by  the  flat  aperture.  The  entire  system  is  in 
a  Sun-Earth  L2  orbit,  which  has  yet  to  be  designed.  The  analysis  of  the 
detector’s  orbit  relative  to  the  aperture  plane  was  beyond  the  scope  of  this 
work. 

At  the  beginning  of  this  task,  we  were  directed  by  the  sponsor  to  inves¬ 
tigate  the  following: 

1.  types  and  fidelity  of  models  used  to  describe  the  relative  orbits  of  the 
seven  spacecraft  forming  the  aperture 

2.  relative  orbits  of  the  seven  spacecraft  forming  the  flat  aperture 

3.  slewing  motions  of  the  flat  aperture 

The  optics  hub  slews  only  in  attitude.  The  entire  system  points  all  over 
the  sky  —  there  is  no  nominal  inertial  orientation.  Also,  we  may  select  any 
distance  between  L2  and  the  optics  hub.  Additionally,  the  question  was  raised 
as  to  whether  the  motions  of  the  telescope  relative  to  the  hub  be  described 
in  cylindrical  coordinates.  This  may  be  easier  to  describe  when  working  with 
the  project  astronomers. 


We  were  unable  to  do  all  of  this  with  the  funding  available  for  this  phase 
of  the  project.  The  work  presented  here  covers  the  depth  needed  for  items  1 
and  2.  We  would  need  to  continue  to  extend  the  development  for  one  tele¬ 
scope  into  that  for  six  telescopes.  A  computer-generated  visualization  of  the 
aperture  plane’s  motion  based  on  our  equations  of  motion  would  be  highly 
useful. 

We  have  concerns  about  being  able  to  arbitrarily  select  the  distance  be¬ 
tween  L2  and  the  optics  hub.  We  suggest  that  the  MAXIM  Pathfinder  mis¬ 
sion  may  best  work  in  a  halo- type  orbit.  A  halo  orbit  is  defined  when  the 
frequency  ratio  of  out-of-plane  to  in-plane  periods  is  a  rational  number.  The 
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three-dimensional  halo  orbit  closes  on  itself  and  is  periodic,  in  contrast  to 
a  Lissajous  orbit  in  which  the  trajectory  does  not  close.  In  particular,  we 
recommend  a  halo  orbit  in  which  the  in-plane  and  out-of-plane  periods  are 
equal. 


In  Section  3  we  form  two  sets  of  differential  equations  of  relative  motion: 
full  nonlinear  and  truncated  nonlinear  equations. 

In  Section  4  we  present  an  analytical  solution  to  these  equations. 

In  Section  5  we  show  simulations. 

In  Section  6  we  summarize  this  report. 


3  Relative  Motion  Differential  Equations 

In  Part  1  of  this  research  [1] ,  the  general  second  order  differential  equations  of 
motion  were  constructed  for  an  object  near  the  Sun-Earth  L2  libration  point, 
using  the  force  model  of  the  classical  circular  restricted  three  body  problem. 
In  this  model,  the  Earth  is  treated  as  being  in  a  circular  orbit  about  the  sun, 
the  spacecraft  mass  is  considered  to  be  negligible  as  compared  to  the  two 
primaries,  and  only  point-mass  gravitational  forces  are  considered. 

For  this  system,  depicted  in  Figure  2,  the  differential  equations  of  motion 
for  an  object  (object  i)  near  the  Sun-Earth  L2  are  given  by 


fl^Xe  +  Di)  H2(xe~D2)\  .  2  . 

- o - t - o -  I  x  +  n  ayx, 


where 


r*  =  vector  from  L2  to  object  i 
Hi  =  solar  Keplerian  constant 

fi2  =  terrestrial  Keplerian  constant  (Earth  +  Moon) 
Pn  =  distance  from  Sun  to  object  i 
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P2i  =  distance  from  Earth-Moon  barycenter  to  object  i 
xe  =  distance  from  system  barycenter  to  L2 
D\  =  distance  from  system  barycenter  to  sun 
D2  =  distance  from  system  barycenter  to  earth-moon  barycenter 
x  =  unit  vector  parallel  to  sun-earth  line  of  syzygy, 
pointing  in  sun-to-earth  direction 
n  =  terrestrial  mean  motion  about  sun  (assumed  constant). 

Let  rh  and  rt  denote,  respectively,  the  vector  from  L2  to  the  hub  and  to 
a  telescope.  Therefore,  if  r  is  the  vector  from  the  hub  to  the  telescope,  the 
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differential  equation  of  motion  for  the  telescope  relative  to  the  hub  is 


r  =  rt  -  rh 
hi 


+ 


P  2 


Pit3  P2t3 


rt  - 


P'1  p2  . 

— o  H - o  )  n/i  + 

Plh  P2h 


hl(xe  +  -Dl)  ,  p2  (^e  D2)  , 

O  I  O  I  X 

Pit  P2t 

Pi(xe  +  Di)  p2(xe-D2) 


—  —p  1 


rt 


rh 


Pit3  Pih3 


p2 


Plh° 

rt  rh 


P2ha 


P2t3  P2h3 


Pi{xe  +  D\) 


pit3  plh3 


X  -  p2(xe  -  D2 ) 


P2t3  P2h3 


(1) 


where  now,  in  general,  the  subscripts  h  and  t  refer  to  the  hub  and  telescope. 


3.1  Series  Expansion  of  Differential  Equations 

The  differential  equations  of  motion  are  now  expanded  in  powers  of  the  dis¬ 
tances  between  L2  and  the  hub,  as  well  as  between  the  hub  and  the  telescope. 
This  is  done  with  the  intention  of  developing  an  analytical  solution  in  terms 
of  these  quantities,  useful  for  performing  a  control  system  analysis. 

First,  the  inverse  powers  of  the  p  magnitudes  are  written  in  terms  of  these 
distances.  If  plh  is  the  vector  from  the  sun  to  the  hub,  then 

Pih  =  (xe  +  D^it  +  rh. 


The  square  of  its  magnitude  is  then  given  by 


Pih"  —  Pih  '  Pih 

=  (xe  +  Di)2  +  rh2  +  2(xe  +  Dx)(rh  ■  x) 
=  (xe  +  DiY  +  T/;2  +  2{xe  +  D\)Xh 

2 


—  ( xe  +  D\) 


1  + 


rh 


Kxe  +  D 1 

where  Xh  is  the  x-component  of  r^.  Then, 


+ 


2  Xh 


Xe  +  D 1 


Plh 6 


(xe  +  D 1] 


1  + 


Th 


xe  +  D 1 


+ 


2  xh 


1  -3/2 


xe  +  D 1 


—  ( Xe  +  D\)  3(1  +  €\h )  3//2, 


56 


where 


A 

£lh  ~ 


Th 


+ 


2  Xh 


. xe  +  D  i  /  a:e  +  -Di 

is  assumed  to  be  less  than  unity.  Using  a  binomial  expansion, 

j_=  1  yH/2\.  * 

Plft3  (xe  +  Di)3  ^  V  k  / 


(xe  +  D  i)3 
In  the  same  fashion,  for  the  telescope, 
1  1 


k= 1  v  7 


Pit3  Oe  +  Ah)3 


where 


A 

Ut  — 


T  t 


k= 1 


+ 


2  .xv 


xe  +  Di)  xe  +  Di 
again  assumed  to  be  less  than  unity. 

Now,  rt  may  be  written  in  terms  of  r/(, .  Using 

rt  =  rh  +  r, 


(2) 


Also, 

Therefore, 

Ut  = 


rt  =  (rft  +  r)  •  (rh  +  r) 
=  rh2  +  r 2  +  2rft  •  r. 

xt  =  xh  +  x. 


r  t 


+ 


2ay 


xe  +  Th  /  a:e  +  D\ 
rh2  +  r2  +  2rh  •  r  2(xh  +  x) 


(xe  +  Dxy 


Xe  +  D  i 


UP 


+ 


2xh 


[(Xe  +  D^2  Xe  +  Dx 
—  Cl  h  +  ^1) 


+ 


r2  +  2r/j  •  r 


+ 


2x 


[(Xe  +  D^2  Xe  +  D  ij 
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where 


^  a  r2  +  2r h  ■  r  2x 

1  (xe  +  Di )2  +  xe  +  L>!  ' 

Again  using  a  binomial  expansion,  the  powers  of  e± t,  as  required  in  Equa¬ 
tion  (2),  are  given  by 

£it.k  —  (ei  h  +  $i)k 


—  Cl  h 


Eh 


s:  £  fc— £ 
>1  Cl/i 


This  expression  may  now  be  used  in  Equation  (2),  giving 


Pit3  (xe  +  Di)3 


(xe  +  Di)3 


i+E 

fc=l 

oo 

1  +  E 


+E  h 


c  €  k—l 
>i  eih 


-3/2V  * 

k  rh 


Eh 


s:  f  fc— l 

h  £\h 


1  1  / — 3/2 

yy  +  (xe  +  o,)3^f  fc 


E  h^M 


Therefore,  for  use  in  Equation  (1), 


Pit3  Pih3  (xe  +  D i)3  ^  V  k 


-3/2 


E  h^M 


Additionally, 


ft  _  _Th_  _  rh  +  r  _ 

Pit3  Pih3  Pit3  Pih3 


Pit3  Pih3  )  +  Pit3' 


An  identical  development  may  be  used  to  form  the  following  relationships 
for  distances  involving  the  earth: 


£ 


P2t3  P'2h3  {Xe  ~  D 2)3  k 


-3/2 


so--" 


(5) 


and 


where 


rt  rh 


P2t3  P2h3 


=  rh 


P2t3  P2h3 


+ 


p2t 


(6) 


A 

^2 h  ~ 


rh 


XP  —  Dn 


+ 


2  Xh 


xe  —  D  9 


and 


a  r2  +  2r h  ■  r  2x 

=  7 - W  xo  + 


(xe~D2)2  xe-D2 

Substituting  Equations  (3)-(6)  into  Equation  (1),  the  telescope  motion 
relative  to  the  hub  is  then  given  by 


r  =  -/i  i 

~  p2 

—  ~  pi 


(rh  +  (xe  +  E>i)x) 
(rh  +  (xe  -  k?2)x) 
rh  +  (xe  +  Di)x 
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Pit3  Pih3J  Pit3 


+ 


P2t3  P2h3 


P2P 


(xe  +  D  i)3 
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k=  1 
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-3/2 

k 


t= i 


£  V 


.  tc  k-t 
>1  tlh 
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(x.  +  Bi)3 
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k=0 


— 3/2 
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e=o 


£ :  w^-1 


P2 


rh  +  (xe  -  D2)± 


xP 


D2)3 
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k= 1 


-3/2 
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i=i 


£  M 


,  £  k-t 
J 2  ^2 h 
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(xe  ~  £2)3 


£ 

k= 0 


-3/2 

k 


i=o 


£  * 


.  te  k-t 
>2  e2h 


(7) 
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3.2  Magnitude  Ordering 

A  magnitude  ordering  system  is  now  employed  in  order  to  truncate  Equa¬ 
tion  (7),  prior  to  forming  the  solution.  It  is  noted  that  the  motion  of  the 
system  takes  place  on  two  separate  distance  scales.  There  is  a  large  distance 
scale,  in  which  the  motion  of  the  hub  about  L2  is  described,  relative  to  the 
motion  of  L2  within  the  context  of  the  three-body  problem.  Then,  there  is  a 
small  distance  scale,  in  which  the  motion  of  the  telescope  about  the  hub  is 
described,  relative  to  the  motion  of  the  hub  about  L2. 


Table  1:  Distance  Ratios  and  Basic  Accelerations 


r/rh 

8.333  x  10~7 

Th/(xe  +  Dx) 

3.971  x  10”3 

rh/(xe  -  D2) 

3.980  x  lO’1 

r  /  (xe  +  Dx) 

3.309  x  10~9 

r/(xe  -  D2 ) 

3.316  x  10”7 

Hx/{xe  +  Dx)2 

5.812  x  10~6  krn/s2 

H2/(xe  -  D2)2 

1.775  x  10~7  krn/s2 

For  ordering  purposes,  consider  the  hub  motion  about  L2  to  be  on  the 
order  of  600,000  km,  and  the  telescope  motion  about  the  hub  to  be  on  the 
order  of  500  m.  Using  the  constants  of  Appendix  A,  the  relative  distances 
are  approximated  by  the  ratios  of  Table  1.  In  the  differential  equations, 
these  ratios  scale  the  basic  acceleration  quantities  which  also  appear  in  the 
table.  Terms  involving  these  basic  accelerations  scaled  by  rh./(xe  +  D\ )  and 
rh./(xe  —  D2 )  are  now  designated  as  being  of  order  1;  terms  involving  the 
basic  accelerations  scaled  by  r/(xe  +  Dx)  and  r/(xe  —  D2)  are  designated  as 
being  of  order  3. 

Retaining  terms  in  Equation  (7)  only  through  order  3,  substantial  algebra 
gives  the  truncated  differential  equations  as 

r  =  A  [— r  +  3xx] 

+  B  [3xrh  +  3xhr  +  (3rh  •  r  -  15xxh)±] 

+  C  [(3rh  •  r  -  13xxh)vh  +  §  (- rh 2  -  hxh2)r 
-  f(2xhrh  •  r  -  7xxh 2  +  xrh2)±] , 
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where 


.  _  hi _ h2 

(xe  +  DXY  (. xe  -  D2f 

d  _  hi  h2 

Oe  +  Di)4  (xe  -  D2)A 

s-<  _  hi _ h2 

(xe  +  Lh)5  (xe  -  Da)5 

Note  that  this  truncation  includes  terms  which  are  linear  in  the  coordinates 
of  r  and  no  more  than  cubic  in  the  coordinates  of  r^.  Terms  involving  A, 
B ,  and  C  are,  respectively,  of  orders  3,  4,  and  5;  lower  order  terms  do  not 
appear. 

The  acceleration  vector  r  may  be  written  relative  to  a  rotating  coordinate 
system  which  rotates  at  the  constant  angular  rate  n  about  the  z-axis  normal 
to  the  ecliptic,  and  with  the  x  direction  as  previously  defined.  This  gives 


(9) 

(10) 

(11) 


r  = 


x  —  2  ny  —  n2x 
y  +  2  nx  —  n2y  , 

z 


where  the  column  vector  notation  is  used  to  indicate  the  xyz  vector  compo¬ 
nents. 


3.3  Linear  Telescope  Motion  about  the  Hub 

Consider  now  the  linear  motion  of  the  telescope  about  the  hub.  From  Equa¬ 
tion  (8),  the  linear  differential  equations  are  given  by 

r  =  A  [— r  +  3a:x] , 


or,  in  component  form, 


x  —  2  ny  —  ( n 2  +  2A)x 
y  +  2  nx  —  (n2  —  A)y 
z  +  Az 


0. 


(12) 


As  would  be  expected,  these  differential  equations  take  exactly  the  same 
form  as  the  linearized  hub  equations  discussed  in  [1],  with  the  same  funda¬ 
mental  frequencies;  the  same  solution  approach  may  be  followed.  Clearly,  the 
out-of-plane  motion  in  the  z  direction  is  decoupled  from  the  motion  in  the 
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xy  plane.  The  z  equation  describes  simple  harmonic  motion,  with  a  solution 
that  may  be  written  as 

z  =  Az  sin [yt  +  ?/?), 

where  v2  =  A. 

Following  the  standard  linear  analysis,  the  in-plane  motion  is  assumed  to 
exhibit  a  solution  of  the  form 

x 

_  y 

Accordingly,  the  fourth-degree  characteristic  equation  for  s  is 

s4  —  (A  —  2  n2)s2  —  ( n 2  +  2  A)  {A  —  n 2)  =  0.  (13) 

Keeping  in  mind  the  relative  magnitudes  of  n  («  0.0172  rad/day)  and  A  (« 
0.0012  rad/day),  the  solution  to  this  equation  has  four  distinct  roots.  One 
root  is  positive  real,  and  corresponds  to  a  divergent  mode: 

-  n 2  +  n2)2  +  (n2  +  2 A) (A  -  n 2). 

A  second,  negative  real,  root  corresponds  to  a  convergent  mode: 

5C  Sd' 

The  remaining  two  roots  are  a  purely  imaginary  conjugate  pair,  correspond¬ 
ing  to  oscillatory  motion  with  natural  frequency  A  given  by 

A2  =  -f  +  n2  +  -n2)2  +  (n2  +  2A)(A-n2). 

Each  mode  shape  is  described  by  the  eigenvector  associated  with  the  corre¬ 
sponding  eigenvalue  s. 

As  in  the  analysis  of  [1],  periodic  motion  occurs  if  the  in-plane  initial 
conditions  are  selected  so  as  to  excite  only  the  oscillatory  linear  modes.  From 
the  solution  of  the  eigenvalue  problem,  placing  this  requirement  upon  the 
initial  conditions  gives 

®(°)  =  jy{  0) 

2/(0)  =  -kAx(0), 
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where 


,  A2  +  n2  +  2  A 
2  An 

Under  this  relationship  among  the  initial  conditions,  the  complete  linear 
solution  for  the  telescope  motion  about  the  hub  may  now  be  written  in  the 
form 


x  =  —  Ax  cos(A t  +  (f) ) 
y  =  kAx  sin(At  +  4>) 
z  =  Az  sin  (yt  +  t/0- 

3.4  Halo  Telescope  Motion 

It  is  desired  that  the  constellation  of  telescope  spacecraft  remain  in  approx¬ 
imately  the  same  specified  plane  over  the  few  days’  duration  required  for 
observations.  To  achieve  this  orientation,  a  halo-type  orbit  of  the  telescopes 
about  the  hub  is  selected.  Such  an  orbit  provides  periodic  motion  in  the 
aperture  plane,  with  the  out-of-ecliptic-plane  fundamental  frequency  equal 
to  A,  the  fundamental  frequency  of  the  xy- planar  motion. 

In  order  to  do  this,  the  inclusion  of  higher-order  forces  is  used  to  ad¬ 
just  the  out-of-plane  fundamental  frequency.  Consider  the  ^-component  of 
Equation  (8): 

A  =  -Az 

+  B(3xzh  +  3  zxh) 

+  C{-Ylxxhzh  +  3 yyhzh  -  6 zxh2  +  § zyh2  +  f zzh2). 

Recall  that  A  =  z/2,  and  define  A  such  that 

A  =  A2  -  z/2. 

Then,  this  differential  equation  may  be  written  as 
A  +  A2z  =  A.S  +  B(3xzh  +  3  zxh) 

+  C(-Y2xxhzh  +  3 yyhzh  -  6 zxh2  +  \zyh  +  | zzh2). 

Here,  the  magnitude  of  A  allows  the  term  A^  to  be  treated  as  a  higher 
order  term,  grouped  with  the  terms  containing  the  coefficient  C.  Using  this 
formulation,  the  linear  solution  now  becomes 

z  =  Az  sin(At  +  'll;), 
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with  the  same  fundamental  period  as  the  in-plane  motion. 

Together  with  the  x  and  y  components  of  Equation  (8),  the  system  dif¬ 
ferential  equations  are 

x  —  2  ny  —  n2x  =  —2  Ax 

+  B{- 6xxh  +  3  yyh  +  3  zzh)  (14a) 

+  C(l2xxh2  -  12 yxhyh  -  12 zxhzh  -  6xyh 2  -  6 xzh2) 
y  +  2  nx  —  n2y  =  —Ay 

+  B(3xyh  +  3  yxh)  (14b) 

+  C(-12xxhyh  +  | yyh2  +  3 zyhzh  -  6yxh 2  +  \yzh2) 
z  +  A  2z  =  +  B(3xzh  +  3zxh) 

/  2  o  2  q  2\ 

+  C  (-12 xxhzh  +  3 yyhzh  -  6 zxh  +  \zyh  +  | zzh  ) . 


4  Lindstedt-Poincare  Development 

The  analytical  solution  to  the  expanded  equations  of  motion  of  Equations  (14) 
is  now  developed  using  a  modified  version  of  the  Lindstedt-Poincare  method. 
This  method  provides  for  the  sequential  solution  of  a  system  of  differential 
equations,  ordered  by  magnitude  of  the  terms,  while  simultaneously  placing 
restrictions  on  the  initial  conditions,  in  order  to  ensure  periodic  motion.  In 
this  development,  the  equations  are  expanded  through  third  order  in  the 
small  quantities,  which  is  defined  as  terms  that  are  at  most  linear  in  the 
motion  of  the  telescope  relative  to  the  hub,  and  quadratic  in  the  motion  of 
the  hub  relative  to  L2. 

Introduce  the  non-linear  frequency  terms  by  series  expansion,  and  change 
the  independent  variable  from  t  to  r,  where 

T  =  COt 


and 

LU  =  1  +  Ui  +  lu2  H - , 

assumed  to  be  an  asymptotic  series,  is  used  to  scale  the  linear  frequency. 
Using  primes  to  denote  differentiation  with  respect  to  r,  the  left  side  of 
Equations  (14)  become 


uj2x"  —  2uny'  —  n2x 
u>2y"  +  2umx'  —  n2y 
oj2z"  +  A2z 
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Recall  that  the  forcing  terms  in  the  differential  Equations  (8)  are  ordered 
such  that  the  terms  involving  A,  B,  and  C  are,  respectively,  of  orders  3,  4, 
and  5.  Now,  for  notational  convenience,  these  terms  are  reordered  as  orders  1, 
2,  and  3.  The  solution  vector  is  assumed  to  take  the  form  of  an  asymptotic 
series,  such  that 


X 

X\  +  x2  +  x3  +  •  • 

y 

= 

2/i  +  2/2  +  2/3  H — 

z 

_  h  1  *2  A  *3  1  ■■ 

where  the  ordering  by  subscript  is  consistent  with  the  reordering  of  the  terms 
of  the  differential  equations;  the  order  of  a  given  term  is  specified  by  the 
subscript. 

This  expansion  is  then  substituted  into  the  differential  equations  of  Equa¬ 
tions  (14),  and  terms  collected  by  order.  The  resulting  first-order  equations 
for  xi,  yi,  and  z\ ,  are  given  by 


x'l  —  2  ny[  —  (n2  +  2  A)x\ 
y"  +  2nx'l  +  (A  —  n2)yi 
z'i  +  A222 


=  0. 


(15) 


Note  that  these  equations  are  identical  in  form  to  the  linear  equations  of 
Equation  (12),  with  A  replacing  v  in  the  z  component  equation. 

The  second  order  equations  contain  contributions  from  the  motion  of  the 
hub  about  L2.  As  shown  by  Richardson  [2],  the  hub  motion  relative  to  L2 
may  also  be  expressed  in  a  series  form: 


Xh 

X\h  +  X2h  +  •  • 

Uh 

= 

yih.  +  y2h.  +  •  • 

.  Zh  . 

Zlh  +  Z2h  +  •  • 

In  the  second  order  equations,  only  x\h,  yih,  and  Z\h  are  included;  at  third 
order,  the  2 h  terms  contribute  as  well. 

The  second  order  equations  for  x2,  y2,  and  z2,  containing  terms  which  are 
linear  in  the  position  of  the  hub,  are  then 


x2  —  2  ny'2  —  (n2  +  2  A)x2 
y2  +  2 nx2  +  (A  -  n2)y2 
z2  +  A2z2 


—2ui(x"  -  ny[)  -  QBxixlh  +  3B(y1ylhh  +  zizlh) 
-2u1(y’(  +  nx[)  +  3B(x1yVl  +  yixlh) 
-2ujiz”  +  35(0)1^1^  +  zixih) 


(16) 
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Here,  the  telescope  motion  terms  on  the  right  side  are  now  assumed  to  be 
known  from  the  solution  to  the  linear  equations. 

Next,  the  third  order  equations  for  x:i ,  y3,  and  z3,  containing  terms  linear 
in  the  2 h  hub  position  terms  and  quadratic  in  the  1  h  terms,  are 


x3  —  2ny'3  —  (■ n 2  +  2  A)x3  = 

—  2u>i(x2  —  ny2 )  —  2tU2(x”  —  ny[)  —  u\x'[ 

+  3B(y2ylh  +  z2zlh  -  2x2xlh  -  2xix2h  +  yxy2h  +  zxz2h) 
+  &C[xi{2x\h  -  y\h  -  z2lh)  -  2ylxlhylh  -  2z1xlhzlh] 
y"±  +  2nx[i  +  (A-  n2)y3  = 

-  2cji (y2  +  nx'2)  -  2u2(vi  +  ra'i)  ~ 

+  3H(x2?/ife  +  y2xih  +  xij/2/»  +  J/i®2/») 


(17a) 


(17b) 


+  3C  [— AxiXihyih  +  l/i  (— 2^:^  +  § y\h  +  +  z\yihZ\h[ 

z'l  +  A2^3  = 

—  2cni^2  —  ( 2uj2  +  uj2)z"  +  Azi  ^ 

+  3B(x2Zih  +  Z2Xlh  +  X!Z2h  +  ZiX2h) 

+  3C  [— 4:X\XihZih  +  yiVihZih  +  z\  (— 2x\h  +  +  |m^)]  • 

Again,  the  right  side  terms  are  assumed  to  be  known  from  the  lower  order 
solutions. 

Note  that  the  left  side  terms  at  each  order  take  an  identical  form.  There¬ 
fore,  a  particular  solution  may  formed  for  a  general  periodic  forcing  function. 
That  solution  may  then  be  applied  to  each  of  the  actual  forcing  terms,  and 
the  complete  solution  given  by  superposition.  Higher  order  homogeneous 
solutions  are  not  required,  as  they  are  identical  in  form  to  the  terms  of  the 
linear  solution,  and  thus  are  already  present  in  the  complete  solution. 
Accordingly,  consider  the  general  case  of  the  forced  system,  given  by 


x"  —  2  ny'  —  ( n 2  +  2  A)x 

*4.cos(gr  +  6) 

y"  +  2  nx'  +  (A  —  n2)y 

= 

B  sin  (qr  +  9) 

z"  +  A2z 

T  sin(pr  +  7) 

This  vector  is  representative  of  all  forcing  terms  which  occur,  where  q  and 
p  are  integers.  Next,  following  the  method  of  undetermined  coefficients,  the 
particular  solution  is  assumed  to  take  the  form 


X 

Rx  cos  (qr  +  9) 

y 

= 

Ry  sin(gr  +  9) 

z 

Rz  sin(pr  +  7) 
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Substituting  into  the  general  system  of  Equation  (18),  and  solving  for  the 
solution  amplitudes, 

2  nqB  —  A(q2  +  n2  —  A) 
qA  +  (A  —  2  n2)q2  —  (A  —  n2)(n 2  +  2  A) 

_  2nqA  - B(q2  +  n 2  +  2 A)  _ 

:y  g4  +  (4  —  2n2)q2  —  (A  —  n2)(n2  +  2A)  1  J 

z  X2  —  p2' 

Note  that  the  denominator  of  Rx  and  Ry  is,  in  fact,  the  characteristic  equa¬ 
tion  of  the  linear  system,  as  seen  in  Equation  (13),  where  s2  has  been  replaced 
by  —  q2.  Therefore,  upon  examination  of  this  particular  solution,  either  q  or  p 
equaling  A  corresponds  to  the  resonance  condition.  In  that  case,  the  forcing 
function  has  the  same  frequency  as  the  homogeneous  solution,  and  so  the 
particular  solution  must  instead  have  an  amplitude  which  is  secular  in  r. 
However,  in  order  for  the  solution  to  be  bounded  for  all  r,  such  secular  terms 
may  not  appear.  Therefore,  restrictions  are  placed  on  the  solution,  such  that 
secular  terms  do  not  arise. 

For  the  case  that  p  —  A,  it  is  required  that  T  =  0.  For  q  =  A,  the 
general  x  and  ^/-component  differential  equations  are  combined  into  a  single 
fourth-order  differential  equation: 

x""  —  {A  —  2  n2)x"  —  {A  —  n2)(n2  +  2  A)x  = 

[2 nBq  +  {A  —  n 2  —  q2)A]  cos (qr  +  6).  (20) 

It  is  then  required  that  the  forcing  amplitude  be  zero,  giving  1 

2  nBq  +  {A  —  n 2  —  q2)A  =  0.  (21) 


4.1  Order  1 


The  solution  to  the  Erst  order  system  of  Equation  (15)  is  once  again  given 
by 


X1 

—Ax  cos(Ar  +  4> ) 

l)\ 

= 

kAx  sin(Ar  +  0) 

.  Zl  . 

Az  sin(Ar  +  0) 

1Of  course,  the  equations  could  instead  be  combined  into  a  y ""  equation,  which  would 
give  a  different  form  of  the  required  relationship.  However,  under  the  condition  that  q 
satisfy  the  characteristic  equation,  the  two  relationships  are  equivalent. 
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where  here  the  solution  is  taken  as  a  function  of  r. 


4.2  Order  2 

Here,  in  the  second  order  equations,  as  previously  mentioned,  x\ h,  Vih,  and 
Z\h  appear.  Again  from  the  earlier  analysis,  the  linear  hub  solution  may  be 
written  as 


Xlh 

- Axh  cos(A ht  +  4>h ) 

yih 

= 

kAxh  sin(A ht  +  (f)h) 

.  z*h  . 

Azh  sin (uht  +  iph) 

Also,  the  hub  frequencies  A h  and  Uh  may  be  written  using  a  nonlinear  asymp¬ 
totic  frequency  scaling  function  u)h,  where 


k>h  —  1  +  U>i  h  +  Cc >2h  +  •  •  •  • 


As  already  seen,  the  linear  hub  frequency  is  the  same  as  that  of  the  telescope; 
therefore  A/,  =  \uh  and  uh  =  vujh.  For  now,  also  assume  that  u>\h  =  cui- 
Typically,  a  Lindstedt-Poincare  analysis  is  employed  for  autonomous  sys¬ 
tems.  Here,  however,  the  forcing  function  is  an  explicit  function  of  time, 
arising  from  the  assumed  solution  for  the  hub  motion.  In  the  linear  hub 
solution  of  Equation  (23),  the  independent  variable  t  must  now  be  written 
in  terms  of  r,  as 


t  =  t  I OJ 

T 

=  r(l  -  oq  -  UJ2  +  cai2  -I - ). 


Then,  in  the  hub  solution, 

A  ht  =  Ac Oht 

=  A(1  +  U>i  +  U>2h  +  '  '  '  )t 

=  A(1  +  UJI  +  u2h  H - )(1  -  LU1  -  lu2  +  uq2  -I - )t 

=  A(1  +  u>2h  -  H - )t. 


Dehne  the  second  order  contribution  to  this  frequency  correction  as 


Au>2  —  U>2h  ~ 


giving 


A  ht  —  A(1  +  Au2  +  •  •  •  )t. 
Similarly,  for  the  out  of  plane  frequency, 

Vht  =  Vh(  1  -  cui  -  +  cui2  H - )t 

=  uh(l  —  - )r, 


where  u)  represents  the  frequency  correction  through  order  2. 

Also,  write  the  hub  phase  angles  <fth  and  ifth  in  terms  of  the  telescope 
angles  as 

<fth  =  (ft  +  A<ft 


and 


ifth  =  'ift  +  A  ift, 


where  A0  and  A^  represent  the  offset  in  phase  angles  between  the  hub  and 
telescope  motion.  Substituting  into  the  linear  hub  solution  of  Equation  (23), 
and  retaining  the  frequency  correction  through  u2  and  u2h, 


Xlh 

—Axh  cos(A(l  +  A u2)t  +  (ft  +  A(ft) 

yih 

= 

kAxh  sin(A(l  +  Auj2)t  +  (ft  +  A(ft) 

.  Zlh  . 

Azh  sin(z//j(l  -  Q)t  +  ift  +  Aift) 

(24) 


This  linear  hub  solution  and  the  first  order  telescope  solution  of  Equa¬ 
tion  (22)  are  substituted  into  the  right  side  of  the  second  order  telescope 
equation,  Equation  (16).  After  resolving  the  angles,  the  three  component 
equations  are  then  written  as 


x2  —  2ny2  —  (n2  +  2 A)x2  =  —  2u;iA(A  —  nk)Ax  cos(Ar  +  (ft) 

-3(1  +  %)BAxAxh 

cos(A(2  H-  -\-  2 (p  + 

-3(1  -ki)BArA  xh  cos(XAu2t  +  Acft) 

+  |  BAzAzh  cos  ((A  -  vh{  1  -  u))t  -  Aift) 

-  lBAzAzh  cos  ((A  +  uh(l  -  u ;))r  +  2ift  +  Aift) 

(25a) 

y2  +  2 nx2  +  (A  —  n2)i/2  =  —  2co»i A( — A: A  +  n)Ax  sin(Ar  +  (ft) 

—  3BkAxAxh  sin(A(2  +  A u2)t  +  2(ft  +  A(ft) 

(25b) 
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z*2  +  A2^  =  2u\X2Az  sin(Ar  +  if;) 

-  ^BAxAzh  sin  ((A  +  z/h(l  -  A))t  +  0  +  ^  +  A^) 
+  |JBAa,Azft.sin((A  -  uh{  1  -  A))t  +  0  -  ^  -  A-0) 

-  |5A2Aa;ftsin(A(2  +  Acu2)t  +  ^  +  4>  +  A0) 

+  ^BAzAxh  sin(AAa;2T  —  tfo  +  </>  +  A0). 

(25c) 

Within  Equations  (25),  each  term  of  the  forcing  function  is  treated  in 
the  manner  of  the  general  approach  given  in  Equations  (18)  and  (19).  First, 
consider  the  potentially  resonant  terms  in  Equations  (25a)  and  (25b).  Using 
the  notation  of  Equation  (18), 

A  =  — 2caiA(A  —  nk)Ax 
B  =  —2ui\(—k\  +  n)Ax 

q  =  A 
e  =  4>. 

The  condition  of  Equation  (21),  to  avoid  resonance  terms,  gives 

— 2uo1XAx[2n{— kX  +  n)  X  +  (A  —  n 2  —  A2)(A  —  nk )]  =  0. 

In  general,  this  condition  can  only  be  satisfied  for  uo\  =  0.  Next,  consider  the 
resonant-type  terms  in  Equation  (25c).  Here, 

T  =  2cuiA2H2 
p  =  X 

7  =  -0. 


The  condition  to  avoid  resonance  is 

2(jU\X~  Az  =  0 

and,  again,  this  can  only  be  satisfied  for  u\  =  0.  Accordingly,  u>i  is  now 
taken  to  be  zero.  Note  that  this  requirement  also  gives  Co  =  002,  which  will 
be  used  from  this  point  forward. 

The  particular  solution  to  Equations  (25)  is  then  built  using  the  method 
of  undetermined  coefficients,  as  previously  described.  The  resulting  solution 
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takes  the  form 


x2  — -  P2i-AxAxh  cos(A(2  +  Au2)t  +  20  +  A0) 

+  p22AxAxh  cos(AAcj2t  +  A0) 

+  p23AzAzh  cos  ((A  -  -  uj2))t  -  A0) 

+  p24AzAzh  cos  ((A  +  z/ft(l  -  cc;2))r  +  20  +  A0) 
y2  =  cr2ij4j.j4j.ft  sin(A(2  +  Auj2)t  +  20  +  A0) 

+  a22AxAxh  sm(XAuj2r  +  A0) 

+  a23AzAzh  sin  ((A  -  i/ft(l  -  cu2))t  -  A0) 

+  (J24AzAzh  sin  ((A  +  z/h(l  -  a ;2))r  +  20  +  A0) 

2;2  =  n21AxAzh  sin  ((A  +  ^(1  -  cu2))r  +  0  +  0  +  A0) 

+  K22A3;A^sin((A  -  z/01  -cn2))r  +  0  -  0  -  A0) 
+  K23AzAxh  sin(A(2  +  Au;2)t  +  0  +  0  +  A0) 

+  K24AzAxh  sin(AAcn2T  -  0  +  0  +  A0), 


(26a) 


(26b) 


(26c) 


where  the  coefficients  are 


3B[—2nXk(2  +  A oj2)  +  (1  +  f  )(A2(2  +  Aw2)2  +  n2  -  A)] 

P'21  ~  A4(2  +  Au2)4  +  (A  -  2n2)A2(2  +  Acn2)2  -  {A  -  n2)(n2  +  2 A) 

3(1  —  y)B(A2Au22  +  n2  —  A) 

P 22  ~  A4Aw24  +  (A  -  2n2)A2Aw22  -  (A  -  n2)(n2  +  2 A) 

_  -§-B((A  -  uh(l  -  uj2))2  +  n2  -  A) 

P 23  (A  -  z/01  -  u2))4  +  (A  -  2n2)(A  -  uh(l  -  cu2))2  -  (A  -  n2)(n2  +  2A) 

_  |.B((A  +  1^(1  —  a;2))2  +  n2  —  A) 

^24  (A  +  z//,(l  -  cu2))4  +  (A  -  2n2)(A  +  uh(l  -  cn2))2  -  (A  -  n2)(n2  +  2A) 


3£[-2nA(2  +  Acn2)(l  +  f)  +  k{  A2(2  +  Acu2)2  +  n2  +  2A)] 

A4(2  +  Acn2)4  +  (A  -  2n2)A2(2  +  Acn2)2  -  (A  -  n2)(n2  +  2A) 

— 6nAAcn2(l  -  h\)B 

A4Acn24  +  (A  —  2n2)A2Au;22  —  (A  —  n2)(n2  +  2A) 

_ 3n(A  -  z^(l  -  u!2))B _ 

(A  -  vh(l  -  (jj2))4  +  (A  -  2n2)(A  -  uh(l  -  cu2))2  -  (A  -  n2)(n2  +  2A) 

— 3n(A  +  vh{  1  -  cn2))i? 

(A  +  z//,(l  -  o;2))4  +  (A  -  2n2)(A  +  ^(1  -  u2))2  -  (A  -  n2)(n2  +  2A) 
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ft  21 

ft22 

ft23 

ft24 


3  B 

2^/i(l  —  CU2)  (2A  +  z//j(l  —  u^)) 

3B 

2^(1  —  o;2)(2A  —  i/h(l  —  u;2)) 

3B 

2A2(1  +  Acu2)(3  +  Aci22) 

3  B 

2A2(1  -  Aw22) 


4.3  Order  3 


The  solution  procedure  at  order  3  is  much  the  same  as  for  order  2,  with 
much  more  complicated  expressions.  The  previously  determined  lower  order 
telescope  motion  is  substituted  into  Equations  (17),  along  with  expressions 
for  the  hub  motion. 

For  now,  the  differential  equations  are  examined  only  for  potentially  res¬ 
onant  terms  —  those  terms  in  which  A  appears  as  the  frequency  in  the  r 
domain.  Analysis  of  these  terms  will  provide  relationships  between  the  ac¬ 
ceptable  values  of  Ax  and  Az,  as  well  as  corresponding  expressions  for  uj2- 
(Further  research  will  allow  the  development  of  the  actual  order-3  solution 
terms.) 

Denoting  the  Cartesian  components  of  these  potentially  resonant  terms 
by  7 Zx,  7 Zy,  and  TZz, 


77.r 

TZy 

7 Zz 


— 2cu2A(A  -  nk)  +  3C  {(k2  -  2 )Axh2  +  Azh 2) 

+  ^B[Axh2{kcj2i  +  ka  22  +  2p2i  +  2^22)  (27c 

+  Azh~(K 21  —  ft22 ) ]  Ax  cos(A r  +  <p) 

2o;2A(A  k  —  n)  +  ^kC\Ax^{3k2  —  4)  +  Az^2} 

+  \BAxh2(-kp2i  +  kp22  -  cr 21  +  O’ 22)  Ax  sin(Ar  +  (p) 

2CU2A"  +  |  C[{k2  —  4  )Axh‘  +  3  Azh~) 

+  |-B(d^/)2(p23  —  P24)  —  Axh2(K23  ~  ft24))  +  A  Az  sin(Ar  +  ip). 

(27c 


(27b) 
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Note  that  7 Zx  and  7 Zy  contain  Ax  as  a  factor,  and  77?  contains  Az. 

Now  examine  the  conditions  to  avoid  resonance,  given  by  Equation  (21), 
and  the  simpler  z-component  condition.  It  is  clear  that  the  coefficients  of 
Equations  (27)  cannot  satisfy  these  conditions  with  any  degree  of  flexibility 
in  the  oscillation  amplitudes  —  either  or  both  of  Ax  and  Az  would  be  required 
to  be  zero.  Therefore,  it  is  necessary  to  pursue  an  alternative. 

It  is  possible  to  increase  the  flexibility  among  the  orbital  parameters  by 
introducing  additional  resonance-inducing  terms.  One  way  to  do  this  is  to 
require  that  the  hub  be  in  a  halo  orbit  about  L2.  This  means  that  Vh  =  Xh, 
which  gives 

Oi(l  —  ^2)  =  A(1  +  Acu2).  (28) 

Making  this  requirement  introduces  the  following  additional  resonant  terms: 


7 70  =  3{C'[cos(Ar  —  A0  +  0  +  A(f>)  —  cos(A  r  +  20  +  A  0  —  0  —  A0)] 

+  ^B[(ka23  +  n23  +  2p23)  cos(Ar  —  A0  +  0  +  A </>) 

+  (/cu24  +  k2a  +  2 p24)  cos(Ar  +  2-0  +  A-0  -  0  -  A0)]} 

Az  AXhAZfl 


77,' 


77? 


(29a) 

=  |{C'/c[sin(Ar  +  2^  +  A-0  —  cj)  —  Acft)  +  sin(Ar  —  A^>  +  0  +  A0)] 

+  2B[(kp23  -  023)  sin(Ar  -  Aip  +  0  +  A 0) 

-  {kp2A  +  a24)  sin (Ar  +  20  +  A0  -  0  -  A0)] } 

AzAXhAZh 

(29b) 

=  | {C[(k2  —  4)  sin(Ar  +  0  +  A0  —  A 0) 

+  ( k 2  +  4)  sin(Ar  —  0  —  A0  +  20  +  A0)] 

+  277 [(p22  -  k2i)  sin(Ar  +  0  +  A0  -  A 0)  (29c) 

-  (P21  +  ^22)  sin(Ar  -  0  -  A0  +  20  +  A0)] } 

AxAX}lAZfl. 


It  can  clearly  be  seen  that  the  addition  of  these  terms  introduces  more  flex¬ 
ibility.  In  particular,  note  that  Ax  appears  in  1ZZ',  while  it  is  not  in  77?; 
similarly,  Az  is  introduced  into  the  x  and  y  direction  terms. 

Next,  in  order  to  allow  these  new  terms  to  combine  with  the  original 
resonance  terms,  it  is  required  that  the  phase  angles  match.  All  in-plane 
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trigonometric  arguments  are  now  required  to  be  A r  +  0;  out-of-plane  argu¬ 
ments  are  required  to  be  A r  +  0.  This  leads  to  the  following  six  relationships: 


cos(Ar  +  20  +  A0  —  0  —  A  0)  =  ±  cos(Ar  +  0)  (30a) 

cos(Ar  —  A0  +  0  +  A0)  =  ±  cos(Ar  +  0)  (30b) 

sin(Ar  +  20  +  A0  —  0  —  A0)  =  ±  sin(Ar  +  0)  (30c) 

sin(Ar  —  A0  +  0  +  A0)  =  ±  sin(Ar  +  0)  (30d) 

sin(Ar  +  0  +  A0  —  A0)  =  ±  sin(Ar  +  0)  (30e) 

sin(Ar  —  0  —  A 0  +  20  +  A  0)  =  ±  sin(Ar  +  0).  (30f) 


From  the  hub  halo  motion  analysis  of  [2],  the  hub  phase  angles  must 
satisfy  the  relationship 

0ft  0ft  2^"’ 

for  arbitrary  integer  j.  With  this  in  mind,  Equations  (30)  lead  to  the  re¬ 
quirement 

0-0=(£+!)tt,  (31) 

for  arbitrary  integer  t. 

Richardson  [2]  also  indicates  that  j  must  be  odd  (1  or  3)  in  order  to  avoid 
enormous  hub  motion  amplitudes  which  would  violate  the  series  expansions 
being  employed.  Therefore,  the  same  assumption  is  made  here,  that  j  must 
be  odd.  Applying  these  relationships  to  Equations  (30)  gives 

cos(Ar  +  20  +  A0  —  0  —  A0)  =  — (— 1)£  cos(Ar  +  0) 
cos(Ar  —  A0  +  0  +  A0)  =  (— 1  )e  cos(Ar  +  0) 
sin(Ar  +  20  +  A0  —  0  —  A0)  =  — (— 1)£  sin(Ar  +  0) 
sin(Ar  —  A0  +  0  +  A 0)  =  (— 1 Y  sin(Ar  +  0) 
sin(Ar  +  0  +  A0  —  A0)  =  (2  —  j)  cos(A r  +  0) 
sin(Ar  —  0  —  A0  +  20  +  A0)  =  —  (2  —  j )  cos(Ar  +  0). 
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Using  these  results  and  augmenting  Equations  (27)  with  Equations  (29), 


TZX 


TZy 


TZz 


— 2u;2A(A  —  nk)Ax  +  3 CAX  (k~AXh2  —  2Axh2  +  A2/,2) 

+  | BAX  [(ka 21  +  k<722  +  2/?2l  +  2 P22)Axh 2  +  (fC2l  —  ^22)A2fr2] 

+  6(— 

+  |(  — l)f5(  — /C(J23  +  ^23  +  2p23  —  k(J 24  ~~  ^24  —  2 p24t)AzAxflAzfl 


cos(Ar  +  4 >) 


(32a) 


2ca2A(A£;  —  n)-4x  +  | kCAx  (3k2  A^2  —  4Ax/i2  +  A2h2) 
+  \BAxAXh2(—kp2i  +  kp22  —  &2 1  +  C22) 

+  |(  — 1)<7-B(^P23  —  ^23  +  kp24  + 


(32b) 


sin(Ar  +  0) 

2ca2A2AL2  +  §CA42  (k2  AXh2  —  4:Axh2  +  3 Az^) 

+  2^Az  [(P23  —  P24:)Azh “  —  (^23  —  K24)Arft2]  +  AA2 
-  6(-l)*Ch4x,WU  (32c) 

+  |(  —  1)^(P22  +  P21  +  ^22  —  ft2l)ArAi:fcAzh 

(2  -  j)(-l)€cos(Ar  +  0). 


The  relationship  between  Uh  and  A,  given  by  Equation  (28)  may  now  be 
used  in  the  p,  cr ,  and  k  functions,  changing  some  of  their  forms  to 


~lB(X2Au22  +  n2  -  A) 

P 23  =  A4Aca24  +  (A-  2n2)A2Aca22  -(A-  n2)(n2  +  2  A) 

|U(A2(2  +  Ata2)2  +  n2  -  A) 

P 24  ~  A4(2  +  Aca2)4  +  (A  -  2n2)A2(2  +  Aca2)2  -(A-  n2)(n2  +  2  A) 

— 3nAAca2-B 

a23  =  A4Aca24  +  (A  -  2n2)A2Aw22  -  (A  -  n2)(n2  +  2A) 

— 3n\(2  +  Auj2)B 

~  A4(2  +  Aca2)4  +  (A  -  2n2)A2(2  +  Aca2)2  -  (A  -  n2)(n2  +  2A) 
35 

^21  ~~  2A2(1  + Aca2)(3  + Aca2) 
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w 

K22  -  2A2(1- Acu22) 

For  use  in  the  equations  at  third  order,  the  p,  a,  and  k  functions  need  only 
be  expanded  through  linear  terms  in  Au;2.  When  linearized,  the  functions 
become 


P21  —  ^-B{[ — 8nkX  +  (2  +  /c2)(4A"  +  n2  —  A)] 

+  [—  |(8A2  +  A  —  2n2)A2[— 8nk\  +  (2  +  k2)(  4A2  +  n 2  —  A)] 
+  4(— nkX  +  (2  +  k2)  A2)]  Au;2} 

—  P21a  +  p2lfeAo;2 
_  3(2  -  k2)B 
p22  ~  2(n2  +  2 A) 

A 

—  P22a 

3  B 

P23  ~  ~2(n2  +  2 A) 

A 

~  P23a 

P24  =  ^-®{(4A2  +  n2  —  A) 

—  |A2[16A4  +  8  (n2  —  A)  A2  +  (A  —  n2)(3n2  +  A)]Au;2} 

=  p24a  +  P24bAcn2 

(72,  =  §B{[-2nA(2  +  k2)  +  /c(4A2  +  n2  +  2A)] 

+  [-|(8A2  +  A  -  2n2)A2[— 2nA(2  +  k2)  +  A;(4A2  +  n2  +  2 A)] 
+  (— nA(2  +  A;2)  +  4/cA2)]  Acn2} 

—  ^la  +  (721{,Ac<;2 


022  — 


023  — 


3nAAca2(2  —  k2)B 
(A  —  n2)(n2  +  2A) 

0226 AcU2 

3nXAu2B 
(A  —  n2)(n2  +  2A) 


—  (7236 Aa;2 


(7 24  =  -f  A5{2  -  i[48A4  +  4(A  -  2n2)A 

—  024a  +  0246 Au;2 


(n2  —  A)(n2  +  2A)]Acn2} 
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ft 21  =  ^7  (3  —  4Au;2) 

—  ft21a  +  ^21bAu2 
35 

^22-^5 

A 

—  ft22  a 

ft23  =  (3  —  4Aa;2) 

—  ft23  a  +  ^23bAu2 

35 

K'24  “  2V 

A 

—  ft23a, 


where 

d  =  16  A4  +  A(A  —  2n2)A2  —  (A  —  n2)(n2  +  2A). 

The  subscripts  6  and  a  respectively  denote  those  terms  involving  and  not 
involving  Auj2. 

Keeping  this  in  mind,  and  recalling  that  Au2  is  defined  as  the  difference 
between  u2h  and  u2)  Equations  (32)  may  be  written  in  the  form 

R-x  =  (c^laAc  +  OiibUJ2Ax  +  OL2aAz  +  Oi2i,U2Az )  COs(Ar  +  0) 

Ry  —  (Pi aAx  +  P\b<jJ2Ax  +  p2aAz  +  p2bUJ2Az)  sin(Ar  +  <fi) 

Rz  —  ('yia.Ax  T  ryn}u)2Ax  T  72a^4z  T  r)2b^2Az) (2  j)(  1)^  cos(At  +  0), 

where  the  a,  P,  and  7  coefficients  are  functions  of  constants  and  the  the  hub 
solution  parameters  u2h,  Axh,  and  Azh.  Using  a  tilde  to  indicate  substitution 
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of  uj-2h  for  Acn2,  these  functions  are 

oil  a  —  3  C  (k~Axh2  —  2  Axh2  +  Azh 2) 

+  | B  [(kcr2i  +  ka-22  +  2p2i  +  2p22)Axh2  +  (ft2 1  ~~  ^22)A^2] 
c*7&  =  — 2A(A  —  nk)  —  |-B  [(fcom  +  /c<7226  +  2p2ib)Axh~  +  fi^iftAl^/,2] 

Qga  =  f(  —  1)£[4C  +  B(  —  k(J2  3  +  ft23  +  2  p23  —  /C(724  —  K24  —  2/?24)] 

«2fe  =  — 1(  —  —  k(J23b  +  «23ft  —  ^246  *"  2p24b)AxhAzh 

fix  a  —  jkC(3k2Axh2  —  4:AXh~  +  Azh2)  +  ^BAXh2(—kp2i  +  &P22  —  <7n  +  ^22) 
(3ib  =  2X(Xk  —  n)  +  \BAxh2{kp2ib  +  cr2i&  —  (x22ft) 

@2a  —  |(  — 1)£-B(^P23  —  0"23  +  ^/?24  +  &24  )AxhAzh 
02 b  =  |(  — l)f-B(°'23b  —  kp24b  ~  ^24 b)AxhAzh 

7l  a  =  (  — 1)^[— 6(7  +  |-B(p22  +  P21  +  ^22  —  ^21)]  AxhAzh 
7l&  2(  3)  B(^K,21b  P21b) Ax}lAzfl 

72a  =  |  C  (k2 Axh2  —  AAxb2  +  3 Azh2) 

+  2Z?  [(p23  —  P24) Az/i.2  —  (k23  —  ^24) A. +  A 

72 b  =  2A2  +  | B  (p24bAzh2  +  K23 bAxh2)  ■ 

(33) 

Recall  that  the  conditions  to  avoid  resonance  are  as  specified  by  Equa¬ 
tion  (21)  and  by  the  requirement  that  the  z  component  of  the  resonant  forcing 
function  have  zero  amplitude.  Here,  these  requirements  now  take  the  form 
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where 


Cia  =  2n/?iaA  +  (A-n-  \2)ala 

(2a  =  2n/32a A  +  (A-n-  X2)a2a  ,  , 

Clfo  =  2n/?ibA  +  (A  -  n  -  X2)alb 

(,2b  —  2n/?2f)A  +  (A  —  n  —  X~)a2b. 

Equations  (35)  and  (36)  provide  two  equations  in  the  three  variables  Ax, 
Az,  and  uj2.  However,  a  third  relationship  may  be  introduced.  Consider  the 
order-1  solution  of  Equation  (22).  Using  the  frequency  and  phase  relation¬ 
ships  of  this  section  gives 


Xi 

— Ax  cos(Ar  +  0) 

yi 

= 

kAx  sin(Ar  +  0) 

.  Zi  . 

_  (2  -  j)(— l)C4zcos(Ar  +  (j))  _ 

with  j  either  1  or  3,  and  (  either  0  or  1.  This  part  of  the  solution  corresponds 
to  the  case  where  the  aperture  plane  maintains  a  roll  angle  (  about  the  y- 
axis.  The  range  of  this  roll  angle  is  related  to  the  integers  j  and  i  in  the 
fashion  shown  in  Table  2. 

Table  2:  Relationship  of  Aperture  Plane  Roll  Angle  to  i  and  j 


J 

e 

C 

1 

0 

(0.1/2) 

1 

i 

(-jt/2,0) 

3 

0 

(-ir/2,0) 

3 

i 

(O.jt/2). 

Therefore,  the  ratio  (2—  j)(— lYAz)/Ax  may  be  taken  as  an  approximation 
to  tan£.  Let  q  represent  Az/Ax.  Assuming  Ax  to  be  nonzero,  Equations  (35) 
and  (36)  may  be  written  in  terms  of  q  as 

7l  a  +  716^2  +  72  aV  +  726^2'h  =  0.  (38) 

and 

Clo  +  Cl  b^2  +  C2  aV  +  (2b^2V  =  0  (39) 
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These  two  equations  may  be  solved  simultaneously  for  lu2  and  77.  The 
combination  of  the  equations  leads  to  a  quadratic,  with  two  sets  of  solutions 
for  uj2  and  77.  Depending  on  the  hub  motion  amplitudes,  the  solutions  for 
77  correspond  to  approximate  aperture  plane  roll  angles  from  —  n/2  to  ir/2. 
The  solution  for  uj2  is  then 


where 


uj2  — 


-B  ±  VB2  -  4 AC 
2 A 


(40) 


A  —  7lbC26  —  72&Clfe 
B  =  7laC26  +  7lbC2a  —  72aClfc  —  72feCla 
C  7l«C2a  72aCltM 


and  the  corresponding  77  is 


V  = 


(la  +  ClfAb  7la  +  7lfe^2 


C2  a  +  (2b^2  72 a  +  726^2 


(41) 


In  turn,  each  value  of  77  is  associated  with  an  infinite  set  of  ( AX,AZ )  pairs, 
each  of  which  is,  in  turn,  associated  with  a  set  of  initial  conditions  for  the 
relative  telescope  motion. 

A  rule  of  thumb  may  be  developed  in  order  to  bound  the  selection  of  Ax 
and  Az.  Examining  the  solution  to  the  linear  telescope  equations,  given  by 
Equation  (22),  it  can  be  seen  that  the  approximate  maximum  distance  of  a 
telescope  from  the  hub  is 


dmax  max  ( \J Ax  T  Az  ,  ks\x ) 
=  Ax  max(y/l  +  r/2,  k). 


(42) 

(43) 


Therefore,  if  the  maximum  distance  is  required  to  be  no  greater  than  a  dis¬ 
tance  dmax,  Ax  may  be  selected  such  that 


Ax  < 


bin 


rnax(-y/l  +  77 2 ,  k )  ’ 


and  then  Az  =  77^3,. 
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4.4  Solution  Summary 

The  solution  is  now  constructed  by  adding  Equations  (22)  and  (26),  and  using 
the  frequency  and  phase  angle  relationships  of  Section  4.3.  This  combined 
solution  through  order  2  is  then 

x  =  —  Ax  cos(Ar  +  0) 


+  ( p-2iAxAXh 

—  P24AzAzh(- 

-1)0 

cos  (A  (2  +  Acu2)t 

+  0/i 

+  0) 

(44a) 

+  (P'22AxAxfl 

+  P2,iAzAZh{- 

-1)0 

cos(AAcn2r  +  0/, 

-0) 

kAx  sin  (Ar  + 

■0) 

+  (cr2i  AxAXh 

—  <j24  AzAZh( 

-1)0 

sin(A(2  +  Ac n2)r 

+  0/i 

+  0) 

(44b) 

+  (022  AxAXh 

—  023  AzAzh( 

-1)0 

sin(AAcn2r  +  4>h  • 

-0) 

[(-1  )eAz  cos 

(Ar  +  0) 

+  (^21  AxAzh  +  K2?,AzAXh(— 1)^)  cos(A(2  +  Acn2)r  +  0^  +  0) 

-  (k22 AxAzh  +  K2iAzAxh(-lY)  cos(AAc n2r  +  <ph  -  0)]  (2  -  j), 

(44c) 

where  r  =  (1  +  cn2f).  The  coefficients  are  as  given  following  Equations  (26), 
with  the  modifications  which  follow  Equations  (32). 

5  Simulation 

5.1  Implementation  Procedure 

Various  representations  of  the  telescope  equations  of  motion  were  coded  in 
MATLAB  to  examine  the  differences  among  the  solutions.  Some  of  the  results 
are  presented  in  this  section.  The  motions  of  the  hub  and  only  one  telescope 
were  simulated  because  of  programming  and  interpretation  simplicity  at  this 
point.  In  the  future,  a  visualization  could  be  developed  based  on  selection  of 
initial  conditions  and  results  from  the  simulations. 

The  selection  of  the  initial  state  is  not  arbitrary  because  the  resulting 
uncontrolled  orbit  path  is  likely  to  be  unsatisfactory.  The  solutions  of  these 
models  show  that  the  hub’s  trajectory  diverges  from  a  closed,  periodic  orbit 
in  typically  150  to  200  days.  Exploration  of  initial  conditions  has  been  a  time 
consuming  portion  of  this  research. 
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Additionally,  the  uncontrolled  telescope  motion  will  diverge  from  the 
nearby  reference  path  about  the  hub. 

Difficulties  in  finding  initial  conditions  that  produced  useful  trajectories 
motivated  short  duration  time  frames.  Assume  that  as  the  aperture  plane 
orbits  the  L-2  point,  a  typical  observation  will  last  only  a  few  days,  before 
reorienting  the  plane  to  observe  some  other  target.  By  simulating  runs  of  20 
days  duration,  small  diverging  motions  from  the  nominal  halo  motion  were 
maintained. 

The  baseline  solution  is  the  numerical  integration  of  the  full  equations 
of  the  hub  subtracted  from  that  of  the  telescope,  as  shown  in  Equation  (1). 
In  essence,  there  are  two  separate  vehicles  orbiting  the  L2  point  with  very 
similar  initial  conditions.  The  slight  differences  in  initial  conditions  are  due 
to  the  placement  of  the  telescope  with  respect  to  the  hub. 

One  cannot  arbitrarily  specify  initial  hub  and  telescope  positions  for  either 
the  full,  or  the  truncated  nonlinear  model,  or  the  analytical  model,  and  obtain 
closed  or  nearly  closed  orbits  about  the  L2  point. 

To  begin  the  process,  implement  the  lengthy  algorithms  shown  in  Sec¬ 
tion  4.3: 

•  begin  with  a  selection  of  Azh,  the  amplitude  of  the  motion  along  the 
2-axis  of  the  hub  with  respect  to  L2 

•  use  a  consistent  set  of  physical  quantities  taken  from  Dunham  and 
Muhonen  [3]  for  the  constants  used  in  computing  A,  B ,  C  of  Equa¬ 
tions  (8)  and  for  A  and  k 

•  compute  the  a,  /3,  7  terms  of  Equations  (33) 

•  compute  the  (  terms  from  Equations  (37) 

•  solve  for  both  of  the  uj2  solutions  from  Equation  (40) 

•  maintain  two  solutions  for  77  (Equation  (41))  based  on  two  uj2  terms 

•  then  Az  =  r]Ax  yields  two  choices 

Table  3  shows  three  combinations  of  amplitude  selection,  selections  of  j, 
and  the  results  of  calculations  based  on  these  selections.  The  application 
of  Equation  (G-l)  from  Richardson  [2]  constrains  the  Axh  corresponding  to 
the  selected  Azh.  (This  equation  is  plotted  in  Richardson’s  Figure  H-l.) 


82 


Table  3:  Summary  of  Inputs  and  Resulting  Values  Used  to  Select  Orbit  Size 
and  Orientation 


Axh 

(km) 

Azh 

(km) 

3 

U->2  + 

V 

r,+ 

r 

(deg) 

(deg) 

Az- 

(m) 

Az  + 
(km) 

280,667 

550,459 

1 

-0.29211 

-0.29185 

1.6115 

1.6147 

-58.18 

-58.23 

80.58 

0.08073 

280,667 

550,459 

3 

-0.29211 

-0.29185 

1.6115 

1.6147 

58.18 

58.23 

80.58 

0.08073 

227,219 

250,000 

1 

-0.26069 

-0.11043 

0.45009 

3.9851 

-24.23 

-75.92 

22.50 

0.1993 

227,219 

250,000 

3 

-0.26069 

-0.11043 

0.45009 

3.9851 

24.23 

75.92 

22.50 

0.1993 

211,126 

1,000 

1 

-0.23166 

-0.07373 

0.00178 

905.98 

-0.1020 

-89.94 

0.08902 

45.30 

211,126 

1,000 

3 

-0.23166 

-0.07373 

0.00178 

905.98 

0.1020 

89.94 

0.08902 

45.30 

The  amplitude  of  the  hub  motion  along  the  z-axis  was  first  chosen  and  then 
the  amplitude  of  the  hub  motion  along  the  x-axis  is  constrained  by  this 
relationship,  which  assures  halo  orbit  solution  via  correct  frequency  selection. 

The  choice  of  Azh  of  550,459  km  is  the  largest  z-value  that  can  be  selected 
and  yield  real  roots  for  u2.  An  intermediate  value  of  250,000  km  and  a  very 
small  value  of  1,000  km  were  also  selected  for  presentation.  The  appendix 
shows  the  results  of  many  calculations  of  this  example. 

There  are  only  two  acceptable  choices  for  j,  which  controls  the  sign  on 
the  ^-position  equation.  The  selection  of  j  =  1  corresponds  to  a  left  roll 
orientation  of  the  aperture  plane;  the  selection  of  j  =  3  corresponds  to  a 
right  roll.  This  tabic  shows  the  results  of  calculations  based  on  fixing  t  =  1 
and  Ax  =  50m.  (Results  for  l  =  0  are  not  shown  because  the  rj  values  would 
be  negative  in  Table  3  and  lead  to  negative  values  of  Az.) 

Because  the  equation  for  u>2  admits  two  roots,  these  two  values  are  carried 
forward  to  see  their  effects  upon  the  calculation  for  rj.  The  superscripts  — 
and  +  refer  to  the  sign  on  the  square  root  in  each  solution.  Two  values  are 
calculated  for  rj  and,  in  turn,  two  values  of  the  arctangent  of  rj  to  yield  the 
approximate  aperture  plane  roll  angle,  £. 

The  resulting  magnitudes  of  the  77  terms,  and  their  corresponding  angle 
values  should  be  highlighted.  As  stated  above,  the  large  value  of  Azh  is 
associated  with  the  case  in  which  the  two  roots  for  uj 2  approach  a  single 
double  root.  Accordingly,  the  resulting  values  of  77,  £,  and  Az  are  also  nearly 
identical. 

For  the  intermediate  value  of  Azh,  the  lo2  values  are  different.  The  roll 
angles,  of  the  two  solutions,  are  quite  different.  Finally,  at  the  smallest 
value  of  Azh,  the  u2  solutions  are  quite  different.  The  roll  angles  of  the  two 
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solutions  are  at  extreme  values  —  based  on  rj~,  the  halo  plane  is  nearly  in 
the  ecliptic  (x-y)  plane  and  based  on  r]+,  the  halo  plane  is  nearly  in  the  y-z 
plane. 

Because  the  user  can  select  either  Azh  or  the  desired  roll  orientation, 
the  user  does  not  have  control  over  the  combination  of  halo  plane  size  and 
orientation.  When  considering  the  roll  orientation  of  the  aperture  plane, 
from  — 7T / 2  to  tt/2,  the  user  should  realize  that  three  months  later,  the  orbit 
will  have  moved  about  the  Sun  such  that  the  initial  roll  angle  appears  as  a 
pitch  rotation.  This  offers  some  orientation  flexibility. 

Again,  the  user  does  not  have  total  control  over  hub  position  and,  of 
course,  the  aperture  plane  location  with  respect  to  L2. 

Continue  with  the  implementation  of  these  algorithms: 

•  compute  the  p,  a,  k  terms  on  pages  76-77 

•  compute  the  analytical  expressions  of  Equations  (44) 

Equations  (44)  are  the  analytical  expressions  for  the  motion  of  the  tele¬ 
scope  with  respect  to  the  hub.  By  taking  the  derivatives  of  these  equations, 
one  can  compute  the  full  state  as  a  function  of  time.  These  equations  may 
be  used  by  trajectory  control  designers. 

Also  at  7  =  0,  Equations  (44)  may  be  used  to  determine  the  initial  state 
for  use  by  two  comparison  solutions  requiring  numerical  integration: 

•  Solution  of  Equation  (1)  -  full  nonlinear  equations  of  motion  of  the 
telescope  with  respect  to  the  hub 

•  Solution  of  Equation  (8)  -  truncated,  to  second  order  in  the  distance  of 
the  hub  from  L2,  nonlinear  equations  of  motion  of  the  telescope  with 
respect  to  the  hub 

5.2  Simulation  to  Compare  Results  of  Models 

Here  is  an  example  that  compares  the  results  of  the  three  solution  methods 
for  calculating  the  motion  of  the  telescope  with  respect  to  the  hub. 

“Full”  refers  to  the  numerical  integration  of  the  separate  equations  of 
motion  for  the  telescope  and  the  hub,  each  spacecraft  with  respect  to  L2,  and 
then  subtracting  the  hub  motion  from  that  of  the  telescope.  This  represents 
the  implementation  of  Equation  (1). 
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Table  4:  Initial  Conditions  of  Telescope  With  Respect  to  Hub  Resulting  From 
Analytical  Solution 


x(0) 

-  64.7796 

m 

</(0) 

0.0 

m 

z(0) 

26.4452 

nr 

i(0) 

0.0 

m/day 

m 

4.4205 

m/day 

m 

0.0 

m/day 

“Truncated”  refers  to  the  numerical  integration  of  the  truncated  equa¬ 
tions  of  motion  for  the  telescope  with  respect  to  the  hub.  This  represents 
the  implementation  of  Equation  (8)  with  the  terms  that  have  a  subscript  h 
referring  to  the  linear  expressions  for  the  motion  of  the  hub.  The  linear  hub 
solution  is  given  by  Equation  (23),  which  was  further  simplified  to  give  a 
halo  orbit  setting  both  frequencies  to  the  same  value  and  setting  both  phase 
angles  to  zero. 

The  initial  conditions  for  both  the  full  and  truncated  numerical  integra¬ 
tions  are  obtained  from  the  calculation,  at  t  =  0,  of  the  analytical  solution. 

“Analytical”  refers  to  the  second  order  solution,  which  is  an  explicit  func¬ 
tion  of  time,  given  by  Equations  (44). 

From  the  intermediate  case  of  Table  3,  select  the  out-of-plane  hub  am¬ 
plitude  from  the  L-2  point  to  be  250,000  km  (Azh).  Using  Equation  (G-l)  of 
Richardson  [2],  the  resulting  value  of  Axh  is  227,219.42  km.  Select  the  row 
with  j  —  1  or  3,  corresponding  to  the  desired  aperture  roll  direction.  The 
roll  direction  is  implemented  in  the  solution  via  Equation  (44c).  The  roll 
angle  direction  is  also  indicated  in  Table  2.  Remembering  that  Table  3  was 
computed  using  £  =  1,  choose  j  —  1  for  a  negative  roll  or  choose  j  =  3  for  a 
positive  roll.  In  this  simulation,  j  =  3. 

Select  the  initial  amplitude  between  the  hub  and  the  telescope  to  be 
Ax  =  50  m,  along  the  rc-axis.  The  Az~  choice  initially  gives  the  telescope  a 
linear  amplitude  of  22.5  m  away  from  the  hub  along  the  z-axis. 

The  analytical  solution  was  first  evaluated  at  t  —  0  to  obtain  the  initial 
conditions  (of  the  telescope  relative  to  the  hub,  Table  4)  for  the  full  and 
truncated  solutions. 

The  initial  conditions  for  the  hub  are  listed  in  Table  5  with  the  equations 
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Table  5:  Initial  Conditions  of  Hub  With  Respect  to  Sun-Earth  L2  Point 


®h(0) 

227,219.42  km 

Vh{  o) 

0.00  km 

zh(  0) 

250,000.00  km 

0 

0.00  rad 

4>h 

0.00  rad 

Xh{  0) 

(X/k)yh(0)  km/day 

Vh{  0) 

— k\xh(0)  km/day 

4(0) 

—  (2  —  j)\Azh  sin(<4)  km/day 

for  x,  y ,  and  z  reproduced  here  for  convenience.  These  were  given  earlier  in 
Section  3.3. 

The  complete  initial  conditions  for  the  telescope  with  respect  to  the  L2 
point  are  given  below: 

®t(0)  =  Xh(0)  +  x(0) 

Vt(  0)  =  Vh(0)  +2/(0) 

*t(0)  =  z„(0)  +  *(0) 

4(0)  =  xh(  0)  +x(0) 

yt(  o)  =  i/h(o)  +  y(o) 

4(0)  =  4(0)  +  i(0) 

The  following  plots  show  the  small  differences  among  the  three  solution 
types  over  a  20-day  duration.  The  plots  show  the  telescope  motion  with  re¬ 
spect  to  the  hub.  The  origin  of  the  coordinate  system  represents  the  position 
of  the  hub.  Scales  are  expanded  for  the  x  and  z  positions. 

Returning  to  the  discussion  that  began  in  the  introduction  to  Section  5.1, 
we  note  that  over  5  days  each  solution  yields  almost  the  same  result.  The 
analytical  solution  is  quite  accurate  for  trajectory  control  analysis.  We  as¬ 
sume  that  as  the  aperture  plane  orbits  the  L2  point,  a  typical  observation 
will  last  only  a  few  days,  before  reorienting  the  plane  to  observe  some  other 
target.  By  simulating  runs  of  20  days  duration,  we  can  see  small  diverging 
motions  from  the  reference  numerical  solution. 

In  Figure  3,  observe  that  both  the  full  and  truncated  solutions  are  essen¬ 
tially  the  same.  Considering  the  scale  of  the  axis  labeled  “x  position”,  the 


- Full - Truncated - Analytical - -Trunc  C=0 


Figure  3:  Compare  3  Solutions  for  Telescope  Motion  Along  the  rc-axis. 


- Full - Truncated - Analytical - -Trunc  C=0 


Figure  4:  Compare  3  Solutions  for  Telescope  Motion  Along  the  y- axis. 
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- Full - Truncated - Analytical - -Trunc  C=0 


Figure  5:  Compare  3  Solutions  for  Telescope  Motion  Along  the  ;j-axis. 


analytical  expression  is  close  to  the  other  solutions  over  5  days.  The  curve 
labeled  “Trunc  C=0”  shows  that  when  the  lowest  order  term  of  the  truncated 
solution  is  set  equal  to  zero,  the  effects  are  not  significant  in  the  x-direction. 

In  Figure  4,  observe  that  all  solutions  are  essentially  the  same  for  the 
scale  shown. 

In  Figure  5,  we  see  that  both  the  full  and  truncated  solutions  are  essen¬ 
tially  the  same.  The  analytical  solution  stays  close  to  these  solutions  over 
5  days.  Remember  the  truncated  solution  has  all  of  terms  A,  B,  and  C  in¬ 
cluded.  The  highlight  is  the  curve  labeled  “Trunc  C=0”  showing  nearly  the 
same  result  as  the  analytical  solution.  Now  the  importance  of  the  “C”  term 
is  evident.  It  indicates  that  a  more  complete  analytic  solution  would  include 
higher  order  terms. 


6  Summary 

This  report  details  the  preliminary  work  describing  the  formation  flying  be¬ 
tween  spacecraft  near  the  Sun-Earth  L2  libration  point,  beginning  with  the 
circular  restricted  three-body  problem  for  the  hub  motion  about  L2. 


The  halo  orbit  was  specified  as  a  desirable  aperture  plane  orbit,  instead 
of  a  Lissajous  orbit,  because  it  is  periodic.  Additionally,  over  a  one  or  two 
day  observation  period,  the  spacecraft  telescopes  of  the  aperture  plane  will 
minimally  separate  due  to  natural  motions. 

The  analytical  solution,  Equations  (44),  for  the  motion  of  a  typical  tele¬ 
scope  relative  to  the  hub  is  presented,  including  terms  which  are  linear  in  the 
hub  motion.  We  anticipate  feedback  from  our  sponsor  concerning  the  utility 
of  these  equations  for  use  in  orbit  control  system  design. 

We  developed  a  full  non-linear  solution  of  hub  motion,  with  respect  to 
L2,  subtracted  from  telescope  motion,  with  respect  to  L2,  as  a  reference 
trajectory. 

We  developed  a  truncated  non-linear  solution  to  the  telescope  motion 
with  respect  to  the  hub.  Both  the  non-linear  model  and  the  truncated  model 
are  solved  by  numerical  integration. 

All  three  models  produce  trajectories  that  diverge  from  a  closed  path  over 
the  course  of  one  orbit  (roughly  200  days).  We  compared  the  three  models 
over  a  shorter  duration  of  20  days  expecting  that  a  science  observation  would 
be  on  the  order  of  a  few  days,  before  reorienting  the  aperture  plane. 

The  analytical  solution  algorithm  is  quite  lengthy.  One  example  was 
presented  and  used  to  compare  the  accuracy  of  the  solution  compared  to  the 
reference  trajectory.  Over  5  days,  the  analytical  solution  is  very  close  to  the 
reference  trajectory. 
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Appendices 

A  Constant  Parameters 


Table  6:  Physical  Constants  [3] 


Hi 

132,712,440,017.987  knvVs2 

R2 

403,503.236  km3/s2 

n 

0.199106385  xlO”6  rad/s 

Tabic  7:  Derived  Constants 


xe 

151,105,099.094445  km 

D\ 

454.84086785372  km 

d2 

149,597,415.850132  km 

Table  8:  Constants  Specific  to  Sun-Earth  L2  Point  [2] 


h 

-14.82882087 

1/(DU2-TU2) 

1 2 

1.673691389 

1/(DU2-TU2) 

-Si 

-0.7444513767 

1/(DU-TU2) 

•S  2 

0.1250471558 

1/(DU-TU2) 

DU 

1.507683382  x  106 

km 

TU 

58.13235527 

days 

Richardson’s  algorithm  ([2],  page  2-31)  for  computing  Axh  given  Azh  and 
for  computing  oo2h  are  given  here: 

A2xh  +  ^A2zh  +  A  =  0 

u2h  =  siA2xh  +  s2A2h 


91 


These  are  the  modifications  used  to  dimensionalize  Richardson’s  algo¬ 
rithms: 


Axh  —  DU  \ 


-£2(Azh/DUY  -  ATU2)/£i 


^2  h  —  s\(AXh/ DU)2  +  s2(Azh/ DU)2 


Where  in  Table  8,  the  constants  “DU”  refers  to  distance  unit  and  “TU” 
refers  to  time  unit. 
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Table  9:  Computed  Values  and  Coefficients 


A 

1.16605228517927  x  10“3 

1 / day2 

B 

5.84853993419721  x  10"10 

l/(km-day2) 

C 

3.86667873919725  x  10"16 

l/(km2-day2) 

A 

3.53850956958284  x  10”2 

rad/day 

k 

3.18712225987377 

z/2 

1.16605228517927  x  10~3 

rad/day2 

A 

8.60527122236636  x  10~5 

1/ day2 

d 

2.56733055729898  x  10~5 

1/ day4 

P'21a 

1.18887100881492  x  10~6 

1  /km 

P21b 

-6.40826960245161  x  10“7 

1  /km 

P22a 

-2.72318375685819  x  10~6 

1  /km 

P22b 

0 

1  /km 

P23a 

-3.33815613931628  x  10”7 

1  /km 

P23b 

0 

1  /km 

P24a 

1.41409729921967  x  10~7 

1  /km 

P24b 

-1.21027861225803  x  10~7 

1  /km 

cr21  a 

6.51772783060505  x  10"7 

1/km 

<021  b 

-7.61518233174130  x  10"7 

1  /km 

<022  a 

0 

1  /km 

<022  b 

-3.81021051605070  x  10"6 

1  /km 

<J23a 

0 

1  /km 

<023 b 

4.67066447286556  x  10"7 

1  /km 

<024 a 

-8.32024709778584  x  10"8 

1  /km 

<024 b 

1.30305530510727  x  10”7 

1  /km 

K-21a 

2.33548302511691  x  10“7 

1/km 

^21  b 

-3.11397736682255  x  10“7 

1/km 

^22a 

7.00644907535074  x  10“7 

1/km 

K-22  b 

0 

1  /km 

^23  a 

2.33548302511691  x  10~7 

1  /km 

K-23  b 

-3.11397736682255  x  10“7 

1  /km 

K-24  a 

7.00644907535074  x  10"7 

1  /km 

^24  b 

0 

1  /km 
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1  Introduction  and  Problem  Definition 


Some  continuing  long-term  goals  of  our  sponsor  are  to  develop  high-fidelity 
equations  of  motion  representing  the  formation  flying  of  a  spacecraft  con¬ 
stellation  near  the  Sun-Earth  L2  point  and  equations  describing  the  relative 
motion  between  these  constellation  members. 

This  is  our  third  report,  which  continues  from  Part  1  [1]  and  Part  2  [2], 
In  this  report,  we  further  the  work  of  Part  2  by  developing  the  elliptical 
restricted  three-body  problem  from  the  previous  work  with  the  circular  re¬ 
stricted  three-body  problem.  Additionally,  this  new  work  includes  the  force 
perturbations  of  lunar  gravitation,  solar  radiation,  and  spacecraft  thrusters. 
The  effects  of  the  elliptical  orbit  of  the  Earth-Moon  system  about  the  sun 
and  the  force  perturbations  are  incorporated  as  additive  perturbations  to  the 
circular  restricted  three-body  problem. 

This  work  builds  from  the  previously  developed  circular  restricted  three- 
body  problem  formulation.  The  familiarity  of  the  formulations  gained  from 
Part  2  is  maintained  here.  We  present  the  development  of  the  full  nonlinear 
baseline  model,  which  includes  these  perturbations: 

•  begin  with  circular  restricted  three-body  problem 

•  then  adding  terms  describing  the  elliptical  orbit  of  Earth  around  the 
sun 

•  then  adding  terms  describing  the  moon’s  motion  about  Earth 

•  then  adding  terms  incorporating  the  force  due  to  solar  radiation 

•  finally  adding  terms  for  spacecraft  thrust 

We  identify  the  magnitude  of  the  contribution  of  each  perturbation  to 
the  solution,  in  order  to  help  determine  when  it  should  be  included  in  the 
computations.  We  identify  substantial  modeling  uncertainties.  In  addition 
to  identifying  the  magnitudes  of  the  effects  of  the  included  perturbations,  we 
identify  the  magnitudes  of  what  is  dropped  in  the  truncation  process. 

The  development  continues  with  the  expanded  form  of  the  full  nonlinear 
baseline  to  sufficiently  high  order  such  as  to  capture  the  relevant  contribu¬ 
tions  to  the  model.  The  presentation  of  the  perturbations  is  organized  in 
similar  fashion  to  the  baseline  presentation  so  the  reader  can  easily  compare 
and  contrast  the  models. 
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This  research  is  loosely  motivated  by  formation  flying  concepts  for  the 
MicroArcsecond  X-ray  Imaging  Mission  (MAXIM)  Pathfinder.  A  concept 
configuration  of  the  MAXIM  mission  is  depicted  and  described  below. 


Figure  1:  Aperture  Formation  of  MAXIM  Pathfinder,  Concept  of  June  2002. 

Seven  spacecraft  shown  in  Figure  1  are  in  the  same  plane  forming  a  flat 
aperture.  The  optical  hub  is  in  the  “center”  of  six  free-flying  spacecraft 
(telescopes).  Dashed  lines  do  not  indicate  a  material  connection,  but  rather 
indicate  random  radial  directions  listed  T\  through  r6.  The  scalar  length  of 
r  ranges  from  100  to  500  meters.  Additionally,  the  six  free-flying  spacecraft 
are  loosely  spaced  60°  from  one  another. 

There  is  also  a  detector  spacecraft  located  approximately  20,000  km  and 
90°  out  of  the  plane  formed  by  the  flat  aperture.  The  entire  system  is  in  a 
Sun-Earth  L2  orbit,  which  has  yet  to  be  designed  by  the  mission  team.  The 
analysis  of  the  detector’s  orbit  relative  to  the  aperture  plane  is  beyond  the 
scope  of  this  work. 
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In  Section  2.1,  we  present  the  circular  restricted  three-body  problem’s 
equations  of  motion,  the  coordinate  frame  that  describes  the  motion,  and 
define  terms  used  in  the  formulation. 

In  each  section  below,  as  appropriate,  we  explain  the  effects  of  those  terms 
retained  and  dropped  from  the  derivation  in  the  course  of  the  truncation 
process. 

•  In  Section  2.2,  we  present  the  elliptical  restricted  three-body  problem, 
which  is  used  as  the  baseline  for  the  subsequent  development. 

•  In  Section  2.3,  we  incorporate  the  effects  of  the  lunar  gravity. 

•  In  Section  2.4,  we  incorporate  the  effects  of  the  solar  radiation  pressure. 

•  In  Section  2.5,  we  incorporate  the  terms  used  to  include  spacecraft 
thrusters. 

•  In  Section  2.6,  we  consolidate  the  equations  into  one  location.  There 
is  further  discussion  to  compare  the  relative  magnitude  of  the  terms  so 
that  they  may  be  turned  on  or  off. 

•  In  Section  3.1,  we  present  the  expanded  form  of  the  circular  full  non¬ 
linear  baseline. 

•  In  Section  3.2,  we  present  the  derivation  of  the  expanded  form  of  the 
elliptical  full  nonlinear  baseline. 

•  In  Section  3.3,  we  incorporate  the  effects  of  the  lunar  gravity. 

•  In  Section  3.4,  we  consolidate  the  equations  into  one  location.  There 
is  further  discussion  to  compare  the  relative  magnitude  of  the  terms  so 
that  they  may  be  turned  on  or  off. 

Because  there  is  no  expansion  of  the  forces  from  solar  radiation  pressure 
and  spacecraft  thrust,  the  contributions  presented  in  Sections  2.4  and  2.5  are 
merely  restated  in  Section  3.4. 

One  important  measurement  of  the  hub-telescope  concept  is  knowledge 
of  the  distance  between  the  spacecraft  to  within  millimeters  or  even  smaller. 
However,  while  many  of  the  equations  of  motion  presented  in  Sections  2  and  3 
have  terms  that  require  numerical  values  which  are  listed  in  readily  available 
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publications,  some  values  of  the  needed  physical  parameters  are  given  to  a 
number  of  significant  figures  with  uncertainties  in  their  knowledge  or  mea¬ 
surement.  A  certain  value  may  be  defined  one  way  in  one  text  and  a  different 
way  in  another  text.  This  will  produce  different  results  in  the  numerical  cal¬ 
culations.  Examples  of  these  uncertainties  applicable  to  the  equations  in  this 
report  are  shown  in  Section  4.  Some  discussion  is  provided  indicating  the 
impact  to  the  simulation  results  due  to  differences  in  the  calculations. 

Section  5  summarizes  the  results  of  Appendix  A,  which  contains  a  detailed 
discussion  of  the  sensitivity  of  the  relative  motion  to  small  errors  in  the 
assumed  position  of  the  hub. 

Appendix  B  fulfills  a  request  by  the  sponsor  to  minimally  address  relative 
motion  at  the  L2  point  of  the  Earth-Moon  system. 

Appendix  C  explains  the  eclipse  geometry  and  corresponding  shadow 
model  used  in  the  application  of  solar  radiation  pressure. 
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2  Equations  of  Motion 

2.1  Circular  Restricted  Three-Body  Problem 

In  Parts  1  and  2  of  this  researdi  [1,  2],  the  general  second  order  differential 
equations  of  motion  were  constructed  for  an  object  near  the  Sun-Earth  L2 
libration  point,  using  the  force  model  of  the  classical  circular  restricted  three- 
body  problem.  In  this  model,  Earth  is  treated  as  being  in  a  circular  orbit 
about  the  sun,  the  spacecraft  mass  is  considered  to  be  negligible  as  compared 
to  the  two  primaries,  and  only  point-mass  gravitational  forces  are  considered. 


t  J 


Figure  2:  Coordinate  Axis  Definition 


For  this  system,  depicted  in  Figure  2,  the  differential  equations  of  motion 
for  an  object  (object  i)  near  the  Sun-Earth  L2  are  given  by 


Hi(xe  +  Di)  p2(xe-D2)\  .  2  . 

- T - t - T7 -  X  +  n  xex, 

Pit*  p2iS  J 
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where 


r*  =  vector  from  L2  to  object  i 
p  1  =  solar  Keplcrian  constant 

/12  =  terrestrial  Keplerian  constant  (Earth  +  Moon) 

Pu  =  distance  from  Sun  to  object  i 
p-2i  =  distance  from  Earth-Moon  barycenter  to  object  i 
xe  =  distance  from  system  barycenter  to  L2 
D 1  =  distance  from  system  barycenter  to  Sun 
D2  =  distance  from  system  barycenter  to  Earth-Moon  barycenter 
x  =  unit  vector  parallel  to  Sun-Earth  line  of  syzygy, 
pointing  in  Sun-to-Earth  direction 
n  =  terrestrial  mean  motion  about  Sun  (assumed  constant). 

The  coordinate  frame  of  Figure  2  is  a  rotating  reference  frame  with  origin 
O  at  the  system  barycenter.  The  rc-axis  points  along  the  Sun-Earth  line  of 
syzygy  and  the  z-axis  is  parallel  to  the  system  angular  momentum;  the  y- axis 
completes  the  dextral  coordinate  system. 

Let  Yh  and  rt  denote,  respectively,  the  vector  from  L2  to  the  hub  and  to 
a  telescope.  Therefore,  if  r  is  the  vector  from  the  hub  to  the  telescope,  the 
differential  equation  of  motion  for  the  telescope  relative  to  the  hub  is 


r  =  rt  -  vh 

Pi 


+ 


p  2 


Pit3  p2t3 
Pl  p2 


Plh3  P2h3 


rt  - 

rh  + 


Pi(xe  +  Di)  p2(xe-D2) 

o  '  o 


Pit- 


P2P 


Pi(xe  +  Di)  |  p2(xe-D2)\_ 

n  T  o  I  X 


—  ~Pl 


r  t 


r  h 


3  3  /  d  - 

pit6  Pih 6 


Pih° 

r  t  r  h 


P2h° 


P2t3  P2h3 


Pi(xe  +  Di) 


Pit3  Pih3 


x  -  p2(xe  -  D2) 


P2t3  P2h3 


(1) 


where  now,  in  general,  the  subscripts  h  and  t  refer  to  the  hub  and  telescope. 
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2.2  Elliptical  Restricted  Three-Body  Problem 

The  extension  to  Equation  (1)  is  now  developed  for  the  case  of  the  elliptical 
restricted  three-body  problem.  First,  it  is  necessary  to  locate  the  point  which 
is  analogous  to  the  libration  point  L2  in  the  circular  restricted  problem.  As 
would  be  expected,  such  a  point  exists;  as  the  Sun-Earth  distance  varies,  the 
location  of  the  point  oscillates  along  the  x-axis. 


2.2.1  Elliptical  Problem  Libration  Point  Analog 

Say  that  there  is  an  elliptical  analog  to  the  L2  point  and  that  its  position 
relative  to  the  (assumed  inertial)  system  barycenter  is  given  by 

Re  =  xe-k  +  ye  y  +  zJl. 

Letting  /  refer  to  the  true  anomaly  of  Earth’s  orbit  about  the  sun,  the 
coordinate  system  has  angular  velocity 

wc  =  A 


which  now  is  considered  to  be  varying  throughout  the  year.  Differentiating 

Re, 


Xe  ~  he 

Xe  ~  he  -  2 he 

~  Pxe  ' 

Re  = 

lie  +  fxe 

,  Re  = 

lie  +  fXe  +  2  fxe 

-  PVe 

Z e 

where  the  column  vector  notation  is  used  to  indicate  the  xyz  vector  compo¬ 
nents. 

The  Newtonian  gravitational  force  per  unit  mass  acting  on  an  object  at 
this  point  is  given  by 


F  e/m  = 


yi(xe  +  Di)  n2{xe-D2 ) 


r  le 


r2e° 


H \Ve  H2Ue 


r  ie° 
H\Ze 


r  le 


where 

rie  =  {xe  +  L>i)x  +  yey  +  zez  and 


r2e° 
HzZe 
T  2e3 


r2e  =  (xe  -  D2)x  +  ye y  +  zez, 
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and  D\  and  D2  now  refer  to  the  time-varying  locations  of  Snn  and  Earth 
along  the  line  of  syzygy. 

Clearly,  as  in  the  circular  restricted  case,  the  z  equation  is  decoupled, 
and  admits  a  solution  ze  =  0.  If,  as  anticipated,  the  desired  solution  involves 
ye  =  0,  the  x  and  y  components  of  the  force-acceleration  equation  become: 


x  :  xe  -  f2xe  + 


hi 


+ 


h  2 


(xe  +  D 1)2  {xe  -  D2)2 


y  :  fxe  +  2 fxe  =  —  ^-\xe2f)  =  0. 
xe  at 


0 


From  the  ^/-component  equation,  x2  f  is  therefore  constant.  However,  con¬ 
servation  of  the  Sun-Earth  two-body  angular  momentum  per  unit  mass  h, 
gives 

D2f  =  h, 

where  D  is  the  varying  Sun-Earth  distance  D\  +  D2.  Of  course,  unlike  the 
case  of  the  circular  restricted  problem,  this  distance  varies  throughout  the 
year.  Substitution  for  /  in  the  ^/-component  equation  gives 

—  constant. 


For  later  convenience,  define  this  constant  in  terms  of  the  constant  7,  such 
that 


Taking  the  positive  root,  this  gives  the  definition  of  7  as 


A 


7  = 


Xe  ~  D2 


D 


In  the  case  of  the  Sun-Earth  system,  7  is  approximately  0.01007824;  xe 
varies  between  1.486  x  108  and  1.536  x  108  km  throughout  the  course  of  the 
year. 


2.2.2  Relative  Motion  Equation  Derivation 

Let  Rj  denote  the  position  vector  from  the  system  barycenter  (point  O)  to 
spacecraft  m,; .  Referring  again  to  Figure  2, 

Rj  =  Xek  +  r ;. 
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As  above,  the  coordinate  system  has  angular  velocity 


UJC  = 


/z- 


Differentiating  R.;, 


R?:  =  xex  +  fxe  y  +  r2 

R;:  =  (xe  -  f2xe)x  +  (2 fxe  +  fxe) y  +  f*. 


The  distance  xe  and  its  time  derivatives  may  be  written  in  terms  of  the 
varying  Sun-Earth  distance  D  and  its  derivatives,  using  the  following  defini¬ 
tions  of  the  constants  7  and  p,  and  the  associated  relationships: 


7  = 
7  +  1  = 


xP  -  Do 


D 

xe  +  D\ 
D 


where 


Then, 


P 


A  Pr2 
P 


1  ~P 


Pi 

P 


p  —  Pl  +  P2- 


D 1  =  pD 
D2  =  (1  ~p)D, 


xe  =  7  D  +  D2 
=  (7  +  1  -  p)D 
Xe  =  (7  +  1  -  P)D 
Xe  =  (7  +  1  -  p)D. 


This  gives  the  acceleration  vector  R,:  as 

R,  =  (7  +  1  -  p){D  -  f2D)x.  +  (7  +  1  -  p)(2fD  +  fD) y  +  r,.  (2) 

Similarly,  two-body  relationships  may  be  used  to  write  D  and  D  in  terms 
of  /  and  /.  Using  the  two-body  equations  of  motion  for  the  sun  (could  use 
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Earth): 


Ri  =  —  Dpk 
=  —pDk 
ujc  =  f  z 

Ri  =  -pDi  -  pfDy 

Ri  =  (- pD  +  pf2D)±  +  (-2 pfD  -  pfD) y 
=  p(f2D^D)x-p(fD  +  2fD)y 
=  -j^x  (two-body  force) 

This  gives  the  component  equations 
X  :  p(/V>  -  £>)  =  || 
y  :  -p(fD  +  2fD)  =  =0, 


which  yield 


D-fD  = 


P  2 


pD 2 


2  fD  +  fD  =  0. 

Substituting  into  Equation  (2),  the  acceleration  becomes 


R*  =  -(7  + 1  -  p) 


p2 


x  +  Ti 


-(7  +  !)^  +  ^]x  +  ^ 
-(7  +  1)^-7^]x  +  rl 


Introducing  the  two-body  forces, 


R; 


-(7  +  !)^ -7^; 


x  +  r,:  = 


pi 


P'2 


■iPu  r>P'2i • 

Pit  P2i 


(3) 


Note  that  the  explicit  appearance  of  the  time  derivatives  of  /  has  now  been 
removed. 
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Next,  the  vectors  pu  and  p2i  may  be  written  as 


Pu  =  (xe  +  Djx  +  rj, 

=  (7  +  1)-Dx  +  Tj 

Substituting  into  Equation  (3), 

R,  = 


and 


P2i  =  (xe  -  D2)±  +  Yi 

=  yDx  +  r  j. 


(7  +  !)^  +  7^ 


x  +  iy 


/il:[(7  +  l)Di  +  iy]  -  -^[7 D±  +  ry] 


Pli 


pi  P'2 


Pli3  P2i3 


Vi  - 


P2i 

/  ,  i  \  hi  ,  P'2 

(7  +  1)—  +7- 


dx. 


Pli°  p2i 

In  the  rotating  frame,  r*  and  its  time  derivatives  are  given  by 


Xi 

u  -  hi 

xi  -  hi  -  2 hi  -  f2Xi 

Vi 

,  r*  = 

yi  +  hi 

,  U  = 

iji  +  f  xi  +  2 f±i  -  /2p* 

.  Z*  . 

Zi 

iy 

Note  that  vectors  pH,  p2i,  and  r.j  are  dependent  upon  the  eccentricity 
of  Earth’s  orbit,  as  is  D.  To  avoid  this  dependence,  redefine  these  vectors 
relative  to  reference  positions  of  Sun,  Earth,  and  L2  along  the  line  of  syzygy. 
Denoting  the  reference  distances  by  an  overbar, 


R;  =  Xek  +  iy 

=  f  Z 

Ri  =  he  y  +  iy 

R;  =  -/27ex  +  hey  +  iy- 


Again  using  the  dehnitions  of  7  and  p,  xe  =  (7  +  1  —  p)D. 
Using  this  form  of  the  acceleration, 


Ri  =  (7  +  1  -  p)(-/2x  +  fy)D  +  Yi 


Pi  .  P2 
— 3  +  — 3  I  r* 
Pu3  P2i3 


(7 +i)4 +U7 

Pli3  P2i3 


dx, 


where 


Pu  =  (7  +  l)£>x  +  Yi  and  p2i='yD±  +  Yi. 


(4) 
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Vectors  r*  and  r  *  take  the  same  forms  as  before,  now  relative  to  the  reference 
L2  location. 

For  the  relative  motion,  let 

r  =  rt  -  rh  =  Rt  -  Kh. 


Then, 


r  =  -fj.  i 


hi  (7  +  1) 


+  h27 


(5) 


as  in  Equation  (1). 


2.3  Lunar  Gravitational  Effects 


This  section  discusses  the  lunar  contribution  to  the  telescope  motion  relative 
to  the  hub.  Terms  corresponding  to  the  gravitational  force  of  the  moon  upon 
the  hub  and  telescope  spacecraft  are  included  in  the  relative  equations  of 
motion  (telescope  relative  to  hub).  The  resulting  contribution  is  then  cast 
as  an  additive  perturbation  to  the  elliptical  restricted  three-body  problem 
equations,  such  that  when  the  lunar  motion  is  ignored,  the  contribution  re¬ 
verts  to  the  lunar  mass  placed  at  the  Earth  position.  The  moon  is  treated 
as  a  point  mass. 

As  shown  in  Figure  3,  let  p3i  denote  the  vector  from  the  moon  to  space¬ 
craft  i,  either  the  hub  or  a  telescope,  and  let  /i3  denote  the  lunar  Keplcrian 
constant.  The  lunar  force  per  unit  mass  on  spacecraft  i  is  then  given  by 


— p3 


(6) 


Therefore,  the  contribution  to  the  equations  of  motion  of  the  telescope  rela¬ 
tive  to  the  hub  is 


Pm  Pih 


KP3t  P3h 

where  the  subscripts  h  and  t  refer  respectively  to  the  hub  and  telescope. 


(7) 
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2.4  Solar  Radiation  Pressure  Effects 

The  work  presented  in  this  section  and  in  Appendix  C  was  prepared  by  Pro¬ 
fessor  David  Richardson  of  the  University  of  Cincinnati,  under  subcontract 
for  this  project  [3]. 

The  pressure  from  solar  radiation  imparts  a  tiny  force  on  a  telescope 
spacecraft.  Depending  upon  spacecraft  design  and  distance  from  the  sun, 
the  force  can  perturb  both  the  spacecraft’s  attitude  and  orbit.  The  force  can 
also  be  harnessed  to  beneficially  propel  the  spacecraft. 

We  present  a  model  to  compute  the  force  on  a  spacecraft,  accounting  for 
the  force  reduction  when  the  spacecraft  orbits  through  the  terrestrial  shadow. 

At  a  distance  p\t  from  the  sun,  the  solar  flux  I  (the  irradiance)  acting  on 
the  spacecraft  is  given  by 

^  Pit 

where 

L  =  3.842  x  1026  watts 
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is  the  solar  luminosity  (total  emitted  radiation). 

If  A  is  the  cross-sectional  area  of  a  spacecraft  of  mass  m  projected  normal 
to  the  spacecraft-Sun  line,  then  the  solar  radiation  force  Fs  per  unit  mass 
acting  on  the  spacecraft  is 

„  CrLA  1.0198  x  lO17^  AT 

bs  =  - 5-  = - ^ -  INI  /unit  mass, 

Anrncpff  mP\t 

where  c  is  the  speed  of  light  and  0  <  Cr  <  2  is  the  parameter  characteristic 
of  the  reflectivity  of  the  spacecraft  surface  facing  the  sun: 

Cr  =  0  translucent, 

Cr  =  1  perfectly  absorbent, 

Cr  =  2  perfectly  reflective. 


For  trajectory  motion  that  passes  through  any  portion  of  Earth’s  shadow, 
the  full  disk  of  the  sun  will  be  partially  obscured.  In  the  vicinity  of  L2, 
this  will  occur  at  distances  normal  to  the  line  of  syzygy  of  approximately 
13,420  km  or  less.  In  such  cases,  the  force  expression  above  must  be  scaled 
by  a  “luminosity  reduction  factor”  o  which  ranges  from  zero  (total  eclipse) 
to  unity  (full  sunlight).  The  appropriate  expression  for  the  force  per  unit 
mass  is  then  written: 


F 


S 


1.0198  x  10  17CrAo- 
mP\t 


Pit 


1.0198  x  10 17CRAa 
m|  K7  +  1  )Tx  +  rt||3 


(7  +  l)Tx  +  rt] .  (8) 


The  calculation  of  the  luminosity  reduction  factor,  cr,  is  presented  in  Ap¬ 
pendix  C. 

Our  solar  model  did  not  consider  the  11-year  solar  cycle,  or  estimate  daily 
variations  of  solar  flux,  or  the  geomagnetic  tail. 


2.5  Spacecraft  Thruster  Effects 

This  section  presents  the  derivation  used  to  incorporate  the  effects  of  body- 
mounted  thrusters  in  our  equations  of  motion. 

First,  we  give  the  body-fixed  acceleration  components  imparted  by  the 
thrusters.  These  terms  are  calculated  as  shown: 

xb  =  Fx/m ,  yb  =  Fy/m,  zb  =  Fz/m, 
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A  F thrust 


Figure  4:  Vehicle  Thrust  in  Body-Fixed  Frame 


where  Fx,  Fy,  Fz  are  the  components  of  the  thrust  Fthrust  in  the  body-fixed 
frame,  as  shown  in  Figure  4;  m  is  the  vehicle  mass. 

Consider  an  arbitrary  alignment  of  a  body-fixed  coordinate  frame  with 
respect  to  the  rotating  x,  y ,  z  frame.  The  body-fixed  Xb ,  yb ,  z&  coordinate 
frame  has  its  origin  at  the  spacecraft’s  mass  center. 

The  components  of  thrust  expressed  in  the  body  frame  must  be  trans¬ 
formed  into  the  rotating  reference  frame  in  order  to  correctly  incorporate 
these  forces  into  the  description  of  the  motion.  One  way  to  express  this 
transformation  is  as  follows: 


X 

xb 

y 

=  T 

Vb 

z 

.  Zb  . 

where  T  is  the  transformation  matrix  formed  as  a  combination  of  Euler  ro¬ 
tations. 

ft  is  somewhat  easier  to  visualize  the  individual  rotations  by  considering 
the  inverse  rotation  description,  rotating  instead  from  the  rotating  reference 
frame  to  that  fixed  in  the  spacecraft: 


Xb 

X 

Vb 

=  T-1 

y 

.  Zb  . 

z 
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Figure  5(a)  shows  the  desired  orientation  of  the  body  frame  with  respect 
to  the  rotating  reference  frame.  However,  by  first  assuming  the  coincident 
alignment  of  both  the  body  frame  and  the  rotating  reference  frame,  we  de¬ 
velop  the  inverse  transformation  matrix  T”1  as  a  combination  of  Enler  rota¬ 
tions  through  a  set  of  Euler  angles: 

1.  first  rotate  angle  +  /)  about  the  Zb  axis  as  shown  in  Figure  5(b), 
where  f  =  f(t-  t0) 

2.  second  rotate  angle  9  about  the  new  orientation  of  the  ijb  axis,  as  shown 
in  Figure  5(c) 

3.  third  rotate  angle  (f>  about  new  orientation  of  the  Xb  axis,  as  shown  in 
Figure  5(d) 

Combine  the  sequence  of  rotations  as  follows  with  a  right-to-left  ordering 
of  the  rotation  matrices: 


'10  O' 

~c(9)  0  -s(9) 

cty  +  f)  s(ip  +  f)  O' 

T~l  = 

0  c(<j>)  s(0) 

0  1  0 

-s{^  +  f)  +  /)  0 

0  -s(<j>)  c(<j>) 

s(6)  0  c{9) 

0  0  1 

=  771(^)T-1((9)771(^  +  /), 


where  the  functions  c  and  s  represent  cosine  and  sine,  respectively;  the  ro¬ 
tation  matrix  T”1  refers  to  the  required  Euler  rotation  matrix  about  the 
4- axis.  Finally,  the  desired  transformation  matrix  T  (from  the  body-fixed 
frame  to  the  rotating  frame)  is  expressed  as  the  inverse  of  T-1: 

T=(T-l)-'=T^  +  f)Ts(e)TM . 

These  expressions  premultiply  the  thrust  force  per  unit  mass  acting  upon  the 
telescope  spacecraft. 

Consider  the  following  example  that  demonstrates  comparable  thrust  and 
solar  forces:  A  500  kg  mass  telescope  spacecraft  with  a  1  mN  thruster  can 
produce  an  acceleration  of  15  km/day2.  By  choosing  a  solar  flux  reflectivity 
parameter  Cm  of  1.5  and  placing  various  surface  areas  normal  to  the  sun,  the 
spacecraft’s  acceleration  due  to  solar  radiation  pressure  is  shown  in  Table  1. 

Contrast  these  low  magnitudes  to  accelerations  due  to  terrestrial  and  so¬ 
lar  gravity  (as  given  by  Equation  5)  of  about  4,800  km/day2.  Of  course,  these 
magnitudes  depend  upon  example.  One  point  is  that  both  the  low  thrust  of 
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(a)  frame  orien-  (b)  ip  +  f  about  Zb 

tat  ions 


(c)  9  about  yb  (d)  <p  about  Xb 

Figure  5:  Frame  Rotation  from  Body  to  Rotating  Frames 


electric  propulsion  and  the  low  pressure  of  solar  radiation  produce  acceler¬ 
ations  very  much  lower  than  the  gravitational  effects.  Another  important 
point  is  that  solar  radiation  pressure  can  be  used  for  orbit  control  because 
the  force  from  solar  radiation  pressure  may  be  comparable  to  that  from  a 
thruster  suite. 

2.6  Summary 

In  this  section,  the  force  equations  derived  in  Sections  2.1  through  2.5  are 
combined.  The  resulting  equation  models  the  elliptical  restricted  three-body 
problem,  incorporating  lunar,  solar  radiation,  and  thrust  perturbations. 
Combining  the  expressions  of  this  section,  given  by  Equations  (4),  (6), 
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Table  1:  Acceleration  clue  to  Solar  Radiation  Pressure  vs.  Surface  Area 


Sun-Facing 

Solar  Radiation  Pressure 

Surface  Area 

Acceleration 

(nr2) 

(km/day2) 

1  («  3  ft  x  3  ft) 

1 

100  («  33  ft  x  33  ft) 

10 

150  («  40  ft  x  40  ft) 

15 

1,000  («100  ft  x  100  ft) 

100 

and  (8),  along  with  the  thrust  expressions  of  Section  2.5,  the  differential 
equations  of  motion  of  the  telescope  are  given  by: 


Rt  =  (7  +  1  -  p)(-/2x  +  fy)D  +  iy 


hi  p2 


Pit* 

P3t 


P2tc 


rt  - 


/  .  p2 

(7  +  1)—  +7- 


P  It* 


P2P 


dx 


—  P3‘  3 

P3P 

1.0198  x  1017CWlcr 


+ 


m  (7  +  l)Px  +  r t 

—  177  -h  i)ux  rtj 

c(ip')  —s(ip')  0 

-  c(9) 

0 

3(9) 

T  0 

0 

s(ip')  clip1)  0 

0 

1 

0 

0  c((p) 

s(<P) 

0  0  1_ 

-s(0) 

0 

c(9)_ 

.0  s((p) 

c(cp)_ 

thrust 


m 


thrust 


where  ip'  —  ip  +  f.  The  corresponding  equations  for  the  hub  are 
R  h  =  (7  +  1  -  p)(-/2x  +  fy)D  +  rh 


pl  P'2 


Pl/i" 

P3I1 


p2ti~ 


r  h 


(  |  n  Pi  |  p2 

(7  +  !)— -7  +7- 


Pi/P 


P2h 


Px 


p3 " 

P3h 

1.0198  x  IO^C'rAct 
m||(7  +  l)Dx  +  iy||3 


(7  +  l)Rx  +  rh]. 


ER3B 

lunar 

SRP 


Combining  the  relative  motion  expressions  of  this  section  (Equations  (5) 
and  (7)),  along  with  the  solar  radiation  pressure  of  Equation  (8)  and  the 
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thrust  expressions  of  Section  2.5,  the  differential  equations  of  relative  motion 
for  a  telescope  spacecraft  are  given  by: 


r  = 


~  p  3 


r« 

Pit 3 
hi  (7 
Pst 


Jh_ 

Pihc 

1) 


p2 


rt 

(ht1 


Plt  PlP 


_r h_ 
P2hC 

+  h27 


p2t  P2h" 


Psh 


,P3t.  PShP 

1.0198  X  10 17CRA(T 
m||(7  +  1  )Dk  +  iy|  |3 
1.0198  x  1017CRAa 
m||(7  +  l)Dit  +  r  ^  1 1 ;3 


[(7  +  l)Dk  +  rt] 
[(7  +  1  )D±  +  rh] 


Ox 


c(V0 

o' 

'  c(0) 

0 

S(ey 

T 

0 

0 

+ 

s(V0 

c(V0 

0 

0 

1 

0 

0 

c{4>) 

-s(4>) 

0 

0 

1 

-s(9) 

0 

c(8)_ 

0 

s{4>) 

thrust 

m 


ER3B 

lunar 

SR  Pt 
SRPh 

thrust 


We  present  the  following  example  to  permit  comparison  of  the  relative 
contribution  of  the  terms  to  the  telescope  motion.  The  initial  conditions  are 
listed  in  Table  2;  they  were  taken  from  the  examples  discussed  in  Part  2  of 
the  research.  As  mentioned  in  that  report,  these  values  were  selected  so  as 
to  excite  only  the  oscillatory  linear  modes. 

Figure  6  presents  the  solution  to  numerical  integration  of  four  different 
force  models  selected  from  the  summary  equation  above.  The  models  were 
applied  to  both  the  hub  and  telescope,  and  the  resulting  state  vectors  were 
differenced  in  order  to  determine  the  position  of  the  telescope  relative  to  the 
hub.  In  each  case,  the  same  force  model  was  applied  to  both  the  hub  and 
telescope.  Additionally,  the  same  value  of  the  reference  distance  D  was  used 
for  all  force  models. 

The  reference  solution  is  represented  in  the  figure  by  the  x-axis.  This 
refers  to  the  circular  restricted  model  case.  The  circular  restricted  model  so¬ 
lution  is  obtained  by  including  the  elliptical  restricted  model  (ER3B)  terms, 
but  with  Earth’s  orbital  eccentricity  treated  as  being  zero;  this  duplicates 
the  model  discussed  in  Part  2.  The  other  models  depicted  in  Figure  6  rep¬ 
resent  the  addition  of  ellipticity  (ER3B),  lunar,  and  solar  radiation  (SRP) 
effects  over  20  days.  It  is  noted  that  the  elliptical  contribution  provides  the 
dominant  perturbation  to  the  circular  restricted  solution.  In  this  example, 
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Table  2:  Example  Initial  Conditions 


hub 

(state  rel. 
to  L2) 

telescope 
(state  rel. 
to  L2) 

telescope 
(state  rel. 
to  hub) 

x(0)  (km) 

-227,219.419 

-227,219.483 

-0.064780 

y{ 0)  (km) 

0.0 

0.0 

0.0 

z(0)  (km) 

-250,000.000 

-249,999.974 

0.026445 

i(0)  (km/day) 

0.0 

0.0 

0.0 

2/(0)  (km/day) 

25,625.039 

25,625.044 

0.004421 

i(0)  (km/day) 

0.0 

0.0 

0.0 

mass  m(kg) 

500 

500 

F thrust  (N) 

0 

0 

sun-facing 
area  A  (m2) 

150 

150 

should  the  viewing  of  a  scientific  target  be  limited  to  no  more  than  10  days, 
the  perturbation  due  to  the  elliptical  contribution  is  less  than  1  m. 

In  the  MAXIM  or  similar  missions,  there  may  or  may  not  be  an  actual 
hub  spacecraft  located  at  the  aperture’s  center.  Regardless,  it  is  necessary 
to  treat  the  hub  as  a  central  reference  point  for  locating  the  positions  of 
the  individual  telescope  spacecraft.  However,  unlike  the  gravitational  forces 
present,  the  solar  radiation  pressure  effect  upon  a  spacecraft  with  actual 
mass  and  area  is  substantial  as  compared  to  that  upon  a  “phantom”  hub. 
Therefore,  for  purposes  of  the  simulation,  the  hub  “spacecraft”  is  treated  as 
having  the  same  physical  characteristics  as  the  telescope.  In  this  manner,  the 
hub  is  maintained  as  an  adequate  reference  for  the  position  of  the  telescope. 

The  effects  of  thrusters  were  not  simulated  here,  due  to  the  vast  uncer¬ 
tainties  of  force  (both  magnitude  and  direction)  and  duration.  We  considered 
trade  studies  involving  thruster  forces  to  be  beyond  the  focus  of  this  investi¬ 
gation. 
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Figure  6:  Effects  of  Perturbations  on  Relative  Distance  —  Full  Equations 
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3  Expanded  Equations  of  Motion 

3.1  Circular  Restricted  Three-Body  Problem 

In  Part  2  [2],  Equation  (1)  is  expanded  through  terms  which  are  linear  in 
the  coordinates  of  r  and  no  more  than  cubic  in  the  coordinates  of  r^.  This 
expansion  takes  the  form 


r  =  A3  [— r  +  3a;x] 

+  A4  [3xrh  +  3xhr  +  (3rh  •  r  -  1 5xxh)±] 

+  d5  [(3rft  •  r  -  15xxh)rh  +  § (- rh 2  -  5xh2)r 
-  y( 2xhrh  ■  r  -  7xxh2  +  xrh2)Z] , 

where  the  constants  are  given  by 

a  _  /b  ,  /f  2 

(xe  +  D1y  +  (xe-D2y 

Terms  involving  Ai  are  considered  to  be  of  order  i;  terms  of  order  lower  than 
3  do  not  appear. 

The  acceleration  vector  r  may  be  written  relative  to  a  rotating  coordinate 
system  which  rotates  at  the  constant  angular  rate  n  about  the  z-axis  normal 
to  the  ecliptic,  and  with  the  x  direction  as  previously  defined.  This  gives 


r  = 


x  —  2  ny  —  n2x 
y  +  2  nx  —  n2y  , 

z 


where  the  column  vector  notation  is  used  to  indicate  the  xyz  vector  compo¬ 
nents. 


3.2  Elliptical  Restricted  Three-Body  Problem 

Consider  the  elliptical  restricted  three-body  equations  of  motion  as  given  in 
Equation  (5).  The  right  side  of  this  set  of  relative  acceleration  equations  may 
be  expanded  as  in  Part  2  [2],  Consider  the  effects  of  various  contributions  to 
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the  magnitude  ordering  scheme.  For  ordering  purposes,  take 


p  =  3.04  x  10~6 
7  =  1.01  x  10~2 

^  =  4.0  x  10~3  (rh  =  600,  000  km) 

L  =  3.2  x  10“9  (r  =  0.5  km). 

A  rough  estimate  may  be  obtained  of  the  contributions  that  the  various 
perturbations  make  to  the  circular  restricted  problem  solution.  Say  that  a 
perturbing  term  may  be  treated  as  modifying  the  linear  frequencies  asso¬ 
ciated  with  the  circular  restricted  problem,  and  consider  the  square  of  the 
perturbing  frequency  to  be  roughly  the  magnitude  of  the  coefficient  of  r  in 
the  perturbing  acceleration.  (Recall  from  Part  1  that  the  linear  periods  in 
and  out  of  the  xy- plane  are  approximately  177.566  days  and  184.002  days, 
respectively.)  Then,  after  90  days,  the  effects  of  the  terms  containing  various 
powers  of  r  and  ry  are  given  in  Table  3,  with  effects  included  of  roughly  20  m 
and  larger. 


Table  3:  Along-Ellipsoid  Effect  of  Sun-Earth  Perturbation  on  Solution  (90 
days) 


perturbing 

term 

effect 

(m) 

re 

16.7 

rrh 

289.3 

rrh2 

118.5 

rr^ 3 

47.4 

4 

rr  h 

18.9 

We  recommend  retaining  terms  which  contribute  down  to  approximately 
20  m  (can  keep  additional  /q  terms  if  convenient).  Retaining  the  rry3  terms 
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results  in  the  expanded/truncated  differential  equations 


r  =  A3  [— r  +  3a:x] 

+  A4  [3.XT/,  +  3xhr  +  (3rh  ■  r  —  ISxxjJx] 

+  A5  [(3rft  •  r  -  15xxh)rh  +  §(rft2  -  5xh2)r 
-  f(2xhrh  ■  r  -  7xxh 2  +  xrh2)± ] 

+  A6  [frh(-xrh2  -  2xhr  ■  rh  +  7xxh2)  +  lr(7xh3  -  3 xhrh2) 
+  ^{.-rh?h  ■  r  +  7xrh2xh  +  7xh2r  ■  rh  -  2lxxh3)it\ . 

Next  consider  the  left  side  of  the  differential  equation: 


"  x  -  fy  2fy  -  f 

2x  1 

r  = 

y  +  fx  +  2fx-  f2y 

Z 

L 

x  —  2  ncy  —  n2x 

J 

"  ~b ~ ~ n^y ~ ~ n2^x ' 

= 

jj  +  2  ncx  -  n2y 

+ 

fx  +  2  (/  -  nc)x  -  (/2  -  n2)y 

z 

0 

(10) 


In  this  expanded  form,  the  first  vector  term  represents  the  acceleration  which 
appears  in  the  circular  restricted  problem  (nc  refers  to  the  mean  motion  of 
the  circular  restricted  Earth  orbit).  The  second  vector  term  gives  the  per¬ 
turbation  which  is  added  by  including  the  elliptic  restricted  effects.  The 
perturbation  term  is  to  be  expanded  in  terms  of  the  eccentricity  e  of  Earth’s 
orbit  (~  0.017).  In  keeping  with  the  earlier  magnitude  ordering,  it  is  esti¬ 
mated  that  it  is  sufficient  to  retain  only  contributions  which  are  linear  in  e, 
because  er  is  approximately  8.5  m. 

Say  that  the  circular  restricted  problem  takes  D  as  being  the  reference 
value  D.  Then,  the  corresponding  mean  motion  is  nc  =  \J fi/D3.  Now,  in 
the  elliptic  restricted  problem,  the  mean  motion  is 


n  Tf  n  =  n  = 

ilc.  ii  i y  -Lymean 


If  D  is  chosen  to  be  the  semi-major  axis  a,  then  n  = 
a(l  +  e2/2),  then 

i,  'V2 

n  =  1  1  +  —  J  nc  «  nc. 
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2 


Expressing  /  and  /  in  terms  of  time, 

•  h  J uafl  —  e2)  na2 a/1  —  e2  /  a 

f  =  -  =  — - - - -  —  -  ~  71  I  - 

■'  D2  D 2  D 2 

where  n  =  \/ /i/a3,  as  given  by  Kepler’s  third 
relationship 

D  =  a(l-e2) 

1  +  e  cos  /  ’ 

■  /I  +  e  cos  f\ 2  . 

f  ~ric[  1  _  g2  )  ~  nc(l  +2ecos/). 

In  terms  of  the  mean  anomaly  Z, 

/  ~  nc(l  +  2e  cos£). 

Differentiating, 

/  ~  —2en2  sin  U 

Using  these  results  in  the  perturbing  vector  of  Equation  (10)  gives 

"  ~(y  ~  2(/  _  -  (/2  -  nc2)a:  ' 

f  x  +  2 (/  -  nc)i  -  (/2  -  nc2)y  = 

0 

—2eyn2  sin  £  —  4eync  cos  £  —  4ermc2  cos  £ 
2exn2  sin  £  +  4einc  cos  £  —  ^eyn2  cos  i 

0 


.ZD 


law.  Using  the  two-body 


3.3  Lunar  Gravitational  Effects 

Consider  the  lunar  force  given  above  by  Equation  (7). 

As  in  the  previous  work,  it  is  desired  that  this  contribution  be  expressed 
in  terms  of  the  position  of  the  telescope  relative  to  the  hub.  Let  the  vector 
r  again  refer  to  the  relative  telescope  position.  Therefore, 

Pit  =  Psh  +  r-  (11) 

As  in  the  earlier  analysis  of  the  effects  of  Sun  and  Earth,  the  contribution  of 
Equation  (7)  may  be  written  as  an  expansion  in  terms  of  r  and  its  compo¬ 
nents.  Accordingly,  a  similar  binomial  expansion  development  is  followed  as 
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that  employed  in  Part  2  of  the  research.  The  square  of  the  magnitude  of  p3t 
is  given  by 


Then, 


P3t 


P3t  '  P3t 

P3h  +  r2  +  2 p3h  ■  r 


P3h 


2p3ft  -  r 

^  9 

P3hl 


1  + 


r  y 

p3h  ) 


+ 


p3t  P3h 

=  P3h  3(1  +  53)  3^2, 


2P3h  ■  r 

P3h2 


-3/2 


where 


Y 

P3hJ 


+ 


2  P 


3  h 


■  r 


P3hz 


is  assumed  to  be  less  than  unity.  Using  a  binomial  expansion, 


Using  this  expansion  in  Equation  (7),  and  using  Equation  (11)  to  substitute 
for  P3t  gives 


Expanding  through  linear  terms  in  r, 

i?  ~  , .  r  i  q  , ,  Psh  r  '  P3h 

Fm  ~  /i3  o  3/i3  o  9 

P3h 6  P3hd  P3h 

ft  is  preferable  to  treat  these  terms  as  an  additive  perturbation  to  the 
equations  of  the  elliptical  restricted  three-body  problem.  Assume  that  the 
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baseline  elliptical  restricted  problem  contains  an  object  with  combined  terres¬ 
trial  and  lunar  mass,  located  at  the  mean  Earth-Moon  barycenter.  Therefore, 
the  expansion  of  Equation  (7)  should  contain  the  lunar  contribution  to  this 
combined  mass  along  with  terms  representing  the  effect  of  the  lunar  motion 
about  the  barycenter.  Accordingly,  Equation  (12)  may  be  rewritten  as 


F  m  ~  ~P  3 


■Hhh  r  •  P-2h  \ 

P2h3  P2h2  ) 


-  /LA 


+  3/i3r  • 


PzhPzh 

P3h5 


p2hP2h\ 

) ' 


(13) 


Because  the  vector  p2h  was  used  in  the  earlier  analysis,  it  is  convenient 
to  write  p3h  as 

Psh  =  P2h  T  hmi 

where  rem  refers  to  the  lunar  position  relative  to  the  mean  barycenter.  The 
square  of  the  magnitude  of  p3h  is  given  by 


This  gives 


P3h  —  P5h  '  P3h 

—  P2h 2  +  r  era2  +  2  p 


2h  '  1  em 


—  P2h 


1+  (  —  )  + 

P2h 


Zp2h  •  re 
P2h2 


r. 


1  _  1 
P3h  P2h 

Using  a  Legendre  polynomial  expansion, 


1+1  —  1  + 

P2h 


Zp2h  '  £e 
P2h2 


-1/2 


-=2-fv^yA(cos5), 

P3h  P2h  “  V  P’2h  ) 


where 


a  A  P2h  '  rem 

COSO  =  - . 

P2h'Y’  em 

Expanding  through  terms  linear  in  rem, 


1  1 

_  _ 

P3h  P2h 


T  em  0 

- COS  O 

p2h 


5 
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giving 


P3h 

1 


Pin 

1 


P3h  Plh 


1  +  3— cos  5 
Plh 

1  +  5 - COS  b 

Plh 


Using  these  representations  in  Equation  (13), 


F 


m 


~p3 


f  r  3p2fer-p2fe\ 

\Plh3  P2h5  J 

3/U  [rP2h  '  I*em  +  r  •  P2hr  em  r  •  remp2h 
~  •  P2hP2hPlh  '  I*  era  /P2h\  /P  2h  • 


(14) 


Examining  Equation  (14),  the  form  of  the  first  term  in  Fm  is  identical 
to  the  corresponding  terrestrial  term  derived  in  the  earlier  work.  The  only 
difference  is  that,  here,  the  mass  coefficient  is  /i3  rather  than  /i2 .  This  term 
represents  the  effect  of  a  lunar  mass  collocated  with  the  terrestrial  mass; 
the  remainder  of  Fm  represents  the  effect  of  the  lunar  motion  about  Earth. 
The  first  term  may  be  added  to  the  earlier  work  simply  by  replacing  ^2 
with  /i2  +  // 3  in  the  relative  equations  of  motion  for  the  elliptical  restricted 
three-body  problem. 

For  the  second  term  of  Equation  (14),  the  Earth-hub  position  vector  p2h 
is  expanded  in  a  fashion  similar  to  that  of  the  earlier  work.  In  the  context  of 
the  elliptical  restricted  problem  using  a  reference  location  of  L2, 


Pih  =  +  rh 


where  7  and  D  are  constants  as  previously  defined.  The  square  of  the  mag¬ 
nitude  of  p2h  is  given  by 


P2h 


P2h  ’  Plh 

(7  Df  +  rh2  +  2y  Dxh 
2 


7  D 


1  + 


/ 

V7  D 


+ 


2 Xh 

7  D 


5 


where  Xh  is  the  x-component  of  r^.  This  gives 

1/P2h  =  (1  +  e/i)1/2/ (7-D)- 
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where 

a  /  rh\2  2xh 
6h~  {'yDj  +7 O' 

Once  again,  powers  of  this  fraction  are  formed  using  binomial  expansions, 
giving 


1  _  1 

P-. ~  (7 DY 


The  negative  powers  of  p2h  that  appear  in  the  second  term  of  Equa¬ 
tion  (14)  are  then  formed  using  the  appropriate  values  of  t  in  this  binomial 
expansion.  Expanding  through  linear  terms  in  77  and  substituting,  Fm  be¬ 
comes 


-p3 


^P2hY  '  P‘2h 


P2h °  P2hu 

3/73  [(  2xxem  T  yy fra  T  zzem')x  T  ( yxem  T  xyem) y 

T  (  ZX fan  \~  X a-  era  )z]/(7£>)4 


-  /i3  [~15xemxhr  +  3rh  •  remr  -  15xremxh 

+  3r  •  rftrem  +  3(r  •  rem)(-5x/lx  +  rh)  +  I05xxemxhit 
-  15xxemrh  -  15xrh  ■  remx  -  15xemr  •  rftx]  /  (7 Df. 


(15) 


Again,  a  rough  estimate  of  the  contributions  that  the  various  perturba¬ 
tions  make  to  the  solution  are  presented.  As  for  the  inclusion  of  the  ellipticity 
in  Section  3.2,  the  perturbing  terms  are  treated  as  modifying  the  linear  fre¬ 
quencies.  After  90  days,  the  effects  of  the  lunar  terms  containing  various 
powers  of  the  relevant  variables  are  given  in  Table  4,  including  contributions 
greater  than  roughly  0.9  m. 


Table  4:  Along-Ellipsoid  Effect  of  Lunar  Perturbation  on  Solution  (90  days) 


perturbing 

effect 

term 

(m) 

era 

2.3 

2 

era 

0.6 

rrhrem 

0.9 
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3.4  Summary 

In  this  section,  the  force  equations  derived  in  Sections  3.1  through  3.3  are 
combined.  The  models  for  solar  radiation  pressure  and  thrust  are  repeated 
here  from  the  equation  of  Section  2.6  for  completeness.  The  resulting  equa¬ 
tion  models  the  expanded  elliptical  restricted  three-body  problem,  incorpo¬ 
rating  lunar,  solar  radiation,  and  thrust  perturbations. 

Combining  the  expressions  of  this  section,  the  expanded  differential  equa¬ 
tions  of  relative  motion  for  a  telescope  spacecraft  are  given  by: 


r  =  A3  [— r  +  3xx\  +  A4  [3xr^  +  3xh r  +  (3rh  ■  r  —  17>xxh)x\ 

+  A5  [(3rft  ■  r  -  15xxh)rh  +  §(rh2  -  5xh2)r 
-  f(2xhrh  ■  r  -  7xxh2  +  xrh2)x] 

+  A6  [f  r h(-xrh2  -  2xhr  ■  rh  +  7xxh 2)  +  \r(7xh 3  -  3 xhrh2) 
+  ^(-rh2rh  ■  r  +  7xrh2xh  +  7xh2r  ■  rh  -  21xxh3)x] 


3 H3  [{-2xxem  +  yyem  +  zzem)x  +  ( yxem  +  xyem) y 
{^ZOCqyti  em  )z]/(7T))4 

l-ki  [ -15xemxhr  +  3rh  ■  remr  -  15xremxh 

+  3r  •  rhrem  +  3(r  •  rem)(-5xhx  +  rh)  +  105xxemxhx 
-  15xxemrh  -  15xrh  ■  remx  -  15xemr  ■  r^x]  /(7-D)5 

1.0198  x  lOl7CRAa  r,  n 

[77  — -TyryT- - — T]3  7  +  1  )Dx  +  rh  +  r 

m||(7  +  l)Dx  +  rh  +  r|p 

1.0198  x  1017Cr71<7  r  ,  .  ,  n 


m||(7  +  l)Dx  +  rh\\3 


(7  +  1)-Dx  +  rh 


cpip')  —  s{ip')  0  c(9)  0  s(9)  10  0 

+  slip')  c(ip')  0  0  1  0  0  c(d)  — s(<£) 

0  0  1  —s(9)  0  c{9)  0  s(</>)  c(</>) 


ER3B 


lunar 


thrust 


where  ip'  —  ip  +  f. 

We  present  the  following  example  to  permit  comparison  of  the  relative 
contribution  of  the  terms  to  the  telescope  motion.  The  initial  conditions  are 
the  same  as  those  presented  earlier  in  Table  2. 

As  in  Figure  6,  Figure  7  presents  the  solution  to  four  different  force  mod¬ 
els  selected  from  the  relative  motion  summary  equation  above.  Where  hub 
position  was  required,  it  was  obtained  from  separate  integration  of  the  full, 
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Figure  7:  Effects  of  Perturbations  on  Relative  Distance  —  Expanded  Equa¬ 
tions 

unexpanded  hub  equations  from  Section  2.6,  using  the  same  force  model  as 
for  the  relative  motion. 

Once  again,  the  reference  solution  is  represented  in  the  figure  by  the  x- 
axis.  This  refers  to  the  circular  restricted  model  case.  The  other  models 
depicted  in  Figure  6  represent  the  addition  of  cllipticity  (ER3B),  lunar,  and 
solar  radiation  (SRP)  effects  over  20  days.  Solar  radiation  is  again  treated  as 
being  applied  to  both  the  telescope  and  to  either  a  real  or  phantom  spacecraft 
at  the  hub,  with  the  same  physical  characteristics  as  the  telescope. 
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4  Modeling  Uncertainties 

4.1  Elliptical  Restricted  Three-Body  Problem 

Very  Similar  Values  for  Different  Definitions.  In  computing  the  value 
of  D,  the  reference  distance  between  the  centers  of  the  sun  and  earth,  we  need 
to  determine  the  mean  distance  between  Earth’s  center  and  the  Earth+Moon 
barycenter. 

The  astronomical  unit  (AU)  is  defined  in  Seidelmann  [4]  as  the  radius  of 
a  circular  orbit  in  which  a  body  of  negligible  mass,  and  free  of  perturbations, 
would  revolve  around  the  sun  in  2rr/k  days,  where  k  is  the  Gaussian  grav¬ 
itational  constant.  This  is  slightly  less  than  the  semi-major  axis  of  Earth’s 
orbit.  Begin  with  the  distance  of  the  AU: 

AU  =  149,  597,  870.000  km  [4,  page  700,  IAU  System] 

=  149,  597,  870.660  km  [4,  page  700,  Best  Estimate] 

Yet  compare  these  values  with  that  from  Dunham  and  Muhonen  [5, 
page  200]: 

AU  =  149,  597,  870.691  km 

The  mean  distance  from  Earth  to  the  sun  also  has  different  values: 

1.0000010178  AU  [4,  page  700,  IAU  System] 

1.00000105726665  AU  [4,  page  700,  Best  Estimate] 

The  calculation  of  the  actual  mean  distance  yields  the  following: 

149,598,022.261  km,  using  the  values  from  the  IAU  System 
149,598,028.825  km,  using  the  values  from  the  Best  Estimate 

Again,  compare  these  values  with  that  from  Dunham  and  Muhonen  [5, 
page  200]: 

Earth+Moon  semi- major  axis  =  1.000001018  AU 

=  149,  598,  023.0  km 

Dunham  and  Muhonen  define  their  value  as  the  mean  distance  of  the 
Earth+Moon  barycenter  from  the  sun.  Contrast  this  value  to  that  of  Seiclel- 
mann,  who  specifies  this  value  as  the  mean  distance  of  only  Earth  from  the 
sun;  however,  the  values  are  very  close  to  one  another. 
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We  continue  using  these  prominent  resources  for  the  eccentricity,  e. 

On  page  700,  Seidelmann  gives  the  value  of  0.016708617  for  the  mean 
eccentricity  of  Earth’s  orbit  about  the  sun.  On  page  200,  Dunham  and 
Muhonen  give  the  value  of  0.01670862  as  the  Earth+Moon’s  eccentricity 
about  the  sun. 

Sensitivity  of  D.  From  the  above,  the  derived  constant  D  depends  on 
two  constants  with  different  values: 

D  =  a(  l  +  e2/2) 

=  149,  618,  904.49  km,  using  IAU  System  for  a  and  Seidelmann  for  e 
=  149,  618,  911.06  km,  using  Best  Estimate  for  a  and  Seidelmann  for  e 
=  149,  618,  905.22  km,  using  Dunham  and  Muhonen  for  both  a  and  e 

These  different  values  were  applied  in  separate  simulations  and  do  not 
substantially  affect  the  analysis  results. 

4.2  Lunar  Gravitational  Effects 

Description  of  the  Lunar  Motion  Around  Earth.  The  moon’s  motion 
around  Earth  is  not  specified  by  an  analytical  model.  The  motion  is  specified 
by  the  software  and  ephcmerides  hies  provided  by  the  JPL.  For  the  needs  of 
this  report,  the  JPL  ephemerides  provide  the  positions  of  Sun,  Earth,  and 
Moon  to  very  high  precision.  The  JPL  ephemerides  are  given  as  blocks  of 
Chebyshev  coefficients,  which,  when  interpolated,  reproduce  the  original  JPL 
numerical  integrations  to  within  1.5  cm. 

The  instructions  for  using  the  JPLEPH.200  ephemeris  hies  state  that  one 
calls  the  JPL  model  with  the  moon  as  the  target  and  Earth+Moon  barycenter 
as  reference  point.  Yet  a  comment  was  found  within  the  FORTRAN  code 
that  the  lunar  state  is  always  geocentric. 
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5  Sensitivity  Summary 

It  is  recognized  that  the  hub  position  knowledge  is  likely  to  be  much  less  pre¬ 
cise  than  that  of  the  relative  motion,  perhaps  on  the  order  of  1  km.  There¬ 
fore,  it  is  useful  to  consider  the  effect  of  this  imprecision  upon  the  telescope’s 
position  relative  to  the  hub. 

Appendix  A  presents  a  detailed  derivation  of  the  variational  equations 
used  to  address  the  sensitivity  of  telescope  motion  to  errors  in  knowledge  of 
the  hub  position. 

The  results  of  the  analysis  clearly  show  that  the  errors  in  telescope  posi¬ 
tion  relative  to  the  hub,  based  on  knowledge  of  hub  position  to  1.7  km  are 
sufficiently  small  that  they  may  be  ignored. 

Several  sample  tests  of  this  behavior  were  conducted.  For  one  test  case, 
the  initial  conditions  of  Table  2  were  used  here  as  a  nominal  set  of  initial 
conditions  for  the  hub  and  telescope. 

The  integration  was  then  performed  with  the  hub  offset  from  its  nominal 
initial  state  by  1  km  in  each  directional  component.  This  was  to  simulate 
an  initial  error  in  hub  position.  As  seen  in  Figure  8,  the  telescope  position 
relative  to  the  hub  differs  from  the  nominal  case  by  millimeters  over  40  days. 

The  simulation  was  also  repeated  with  the  initial  hub  position  offset  from 
the  nominal  state  by  17.3  km  (10  km  in  each  direction).  Figure  9  demon¬ 
strates  the  roughly  ten-fold  increase  in  error  over  the  previous  example. 
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Delta  X - Delta  Y - Delta  Z - Mag 


Figure  8:  Relative  Motion  Error  Caused  by  1.7  km  Initial  Hub  Position  Error 

6  Summary  and  Conclusion 

This  report  details  our  further  work  describing  the  formation  flying  between 
spacecraft  near  the  Sun- Earth  L2  libration  point,  beginning  with  the  circular 
restricted  three-body  problem  for  the  hub  motion  about  L2. 

Continuing  from  our  previous  works,  these  analyses  develop  the  elliptical 
restricted  three-body  problem  from  previous  work  with  circular  problem.  We 
built  on  our  familiarity  with  the  circular  problem  to  address  the  following  as 
perturbations  upon  the  circular  problem: 

•  elliptical  orbit  of  Earth- Moon  about  Sun 

•  lunar  gravitational  effects 

•  solar  radiation  pressure  effects 

•  thrusters  on  vehicle 

These  were  incorporated  as  additive  perturbations  to  the  circular  re¬ 
stricted  three-body  problem  with  expansions  of  varying  levels  of  fidelity.  We 
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Figure  9:  Relative  Motion  Error  Caused  by  17.3  km  Initial  Hub  Position 
Error 

develop  two  forms  of  models  of  the  full  derivation: 

•  full  nonlinear  baseline 

•  expanded  form  of  full  nonlinear  baseline  to  appropriate  order 

Section  2  presented  the  derivation  of  the  full  nonlinear  baseline  with  the 
perturbations.  The  equations  were  implemented  through  their  coding  in  a 
MATLAB  simulation.  One  example  was  presented  using  the  discovery  from 
Part  2  of  valid  initial  conditions  that  excite  only  oscillatory  motion  in  the 
linear  modes.  The  results  shown  in  Figure  6  demonstrate  the  dominant 
perturbation  due  to  the  elliptical  motion  of  Earth  about  the  sun.  In  this 
example,  the  perturbations  due  to  lunar  and  solar  radiation  pressure  are 
negligible.  Overall,  the  figure  indicates  that  during  the  typical  5-day  science 
observation  of  this  example,  these  perturbations  contribute  less  than  1  m 
of  relative  position  difference.  The  proposed  electric  propulsion  should  have 
little  trouble  maintaining  position. 
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Section  3  presented  the  expanded  circular  and  elliptical  portions  of  the 
full  model  through  terms  which  arc  linear  in  the  coordinates  of  the  telescope 
position  relative  to  the  hub  and  no  more  than  cubic  in  the  coordinates  of 
the  hub  position  relative  to  L2.  The  derivation  describes  the  magnitude  of 
the  various  terms  and  why  some  were  truncated.  Additionally,  the  lunar 
gravity  model  was  expanded  and  some  terms  truncated  due  to  very  small 
perturbations  upon  the  motion.  As  before,  the  perturbations  due  to  the 
elliptical  motion  of  the  Earth  about  the  Sun  are  dominant.  Again,  lunar  and 
solar  radiation  pressure  effects  are  negligible. 

Modeling  uncertainties  were  discussed  in  Section  4.  One  type  of  uncer¬ 
tainty  is  that  the  same  listed  number  can  be  found  to  have  different  defini¬ 
tions  between  two  popular  references.  Calculations,  of  course,  give  slightly 
different  results;  however,  they  do  not  substantially  affect  the  results. 

Section  5  explains  the  results  of  a  longer  derivation  given  in  Appendix  A. 
During  the  18-months  of  part-time  work  on  this  report,  our  sponsor  requested 
that  we  investigate  the  sensitivity  of  telescope  motion  to  errors  in  hub  posi¬ 
tion.  The  analysis  and  simulation  show  that  errors  in  the  knowledge  of  hub 
position  of  1  km  and  10  km  (in  each  of  the  three  position  components)  yield 
millimeters  of  error  in  telescope  position. 

Appendix  B  presents  a  back-of-the-envelope  look  at  the  relative  motion 
between  a  hub  and  telescope  at  the  Earth-Moon  L2  point.  Our  sponsor  made 
a  comment  that  future  work  may  be  redirected  to  this  vicinity.  Our  work 
addressed  the  effects  of  the  largest  contributions  to  the  elliptical  restricted 
problem.  The  perturbation  caused  by  solar  gravitation  was  found  to  be 
negligible. 

Appendix  C  presents  a  derivation  of  the  luminosity  reduction  factor  in 
the  context  of  the  solar  radiation  pressure  model.  The  luminosity  reduction 
factor  accounts  for  shadowing  effects  due  to  a  partial  eclipse. 


In  conclusion,  based  on  earlier  literature  searches,  we  believe  this  new 
work  is  unique  because  it  describes  the  primary  perturbations  to  the  de¬ 
scription  of  relative  motion  between  nearby  spacecraft.  We  verified  that  the 
effects  of  the  elliptical  motion  of  the  Earth  about  the  Sun  is  the  dominant 
perturbation  to  the  circular  restricted  three  body  problem.  Contributions 
due  to  lunar  gravity  and  solar  radiation  pressure  are  nearly  negligible  in  the 
chosen  examples. 
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Appendices 

A  Sensitivity  of  Telescope  Motion  to  Errors 
in  Hub  Position 

This  appendix  demonstrates  the  low  sensitivity  to  errors  in  hub  position,  of 
the  telescope  motion  relative  to  the  hub.  The  effects  of  1  km  hub  position 
errors  upon  the  relative  telescope  motion  are  presented.  The  model  used  is 
the  elliptical  restricted  model  as  presented  in  Section  3.2. 

A.l  Derivation  of  Variational  Equations 

Consider  the  second-order  differential  equation  of  motion  of  a  telescope  rel¬ 
ative  to  the  hub,  assumed  to  take  the  form  of  the  earlier  expanded  and 
truncated  form: 

r  =  A3  [— r  +  3a:x] 

+  A4  [3xrh  +  3a+r  +  (3rh  ■  r  —  15 xxh)x] 

+  A5  [(3rA  •  r  -  15xxh)rh  +  § (rh2  -  5xh2)r 
-  f(2xhrh  ■  r  -  7xxh2  +  xrh2)x ] 

+  A  [  frh(-xrh 2  -  2xhr  ■  rh  +  7xxh2)  +  lr(7xh3  -  3 xhrh2) 

+  ^(-rh^h  ■  r  +  7xrh2xh  +  7xh2r  ■  rh  -  2lxxh3)x] , 


where 

r  =  [x  y  z ]T 


x  —  2  ncy  —  n2x 

2 eyn2  sin  l  —  4eync  cos  £  —  4 exn2  cos  i 

r  = 

y  +  2 ncx  -  n2y 

+ 

—2 exn2  sin  £  +  4ei;nc  cos  £  —  4 eyn  2  cos  £ 

z 

0 

hi  y-2 

(1  +  7  )nDn  7  nDn' 
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and 


r  =  telescope  position  relative  to  the  hub 
rh  =  hub  position  relative  to  L2 
x,y,z  —  Cartesian  components  of  r 
y i  =  solar  Keplerian  constant 

/i2  =  terrestrial  Keplerian  constant  (Earth  +  Moon) 

D  =  reference  Sun-Earth  distance 
xe  =  distance  from  system  barycenter  to  L2 
D\  =  distance  from  system  barycenter  to  Sun 
D2  =  distance  from  system  barycenter  to  Earth-Moon  barycenter 
7  —  (xe  —  D2)/D 

nc  =  circular  restricted  terrestrial  mean  motion 
t  =  mean  anomaly  of  Earth  orbit. 

The  rotating  coordinate  system  is  defined  to  have  rc-axis  along  the  Sun-Earth 
line  of  syzygy  and  z-axis  normal  to  Earth’s  orbit  about  the  sun. 

The  vector  differential  equation  may  be  written  in  the  form 

fx 

=  f  =  fy  ■ 

.  /*  . 

The  vector  function  f  represents  the  acceleration  relative  to  the  rotating 
frame.  This  function  is  linear  in  x,  y,  z,  and  their  derivatives. 

Say  that  there  is  a  nominal  solution  r„  associated  with  a  nominal  hub 
position  rhn.  It  is  desired  to  estimate  the  sensitivity  of  the  relative  telescope 
motion  to  errors  in  rh}  which  are  on  the  order  of  1  km. 

Define  the  state  vector  x  as 

x  =  [  x  y  z  x  y  z  ]T . 

The  time  derivative  of  this  state  is  then 

x  =  [  x  y  z  fx  fy  fz]T. 

Referring  to  this  vector  as  the  function  F(x,  r^),  the  differential  equation 
now  takes  the  first-order  form 

x  =  F.  (16) 
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where  the  nominal  solution  satisfies 


=  F, 


(17) 


The  solution  x  is  now  considered  to  be  xra  as  perturbed  by  an  error  in 
Yh-  If  is  desired  to  examine  the  effect  of  this  error  upon  x.  Writing  x  as 
a  linearized  Maclaurin  series  in  this  error  (or,  equivalently,  a  Taylor  series 
about  the  nominal  hub  trajectory)  gives 


x  =  x,; 


where 


<9x 

drh 


<9x 

drh 


(rh  -  rhn) , 


(18) 


dx 

dx 

dx 

dxh 

dyh 

dzh 

dy 

dy 

dy 

dxh 

dyh 

dzh 

dz 

dz 

dz 

dxh 

dyh 

dzh 

dx 

dx 

dx 

dxh 

dyh 

dzh 

dy 

dy 

dy 

dxh 

dyh 

dzh 

dz 

dz 

dz 

dxh 

dyh 

dzh 

=  U. 


The  subscript  n  refers  to  evaluation  using  the  nominal  trajectories;  Xh,  Uh, 
and  Zh  are  the  Cartesian  components  of  the  hub  motion  relative  to  L2. 

It  is  the  quantities  of  the  elements  of  this  matrix  U  that  are  of  interest, 
in  particular,  the  partial  derivatives  of  the  position  components.  Once  these 
elements  are  determined,  the  effect  of  the  hub  position  error  upon  the  relative 
telescope  motion  may  be  approximated  by 


dx  dx  dx 
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Because  the  components  of  the  hub  position  error  are  presumed  to  be  no 
greater  than  1  km,  it  is  sufficient  to  examine  the  position  partials  within 
U ;  they  are  treated  as  being  scaled  by  unity  to  form  the  relative  position 
variations. 

A  Maclaurin  series  may  be  formed  as  well  for  the  derivative  of  the  state, 


as 


x  =  x»- 


<9x 

drh 


(rh  -  rhn ) . 


This  derivative  may  instead  be  formed  by  differentiation  of  Equation  (18), 
giving 


•  (rh  -  rhn)  + 


<9x 

drh 


jt  {th  ~ rta)  - 


(19) 


It  is  assumed  that  the  change  in  hub  position  error  is  of  higher  order  than 
the  error  itself.  Therefore,  these  two  equations  indicate  that 


dx 

d  I 

(  dx 

drh 

n  ~  dt  ' 

\  drh 

(20) 


The  function  F  may  also  be  written  as  a  linearized  series,  giving 


F 


Fn  + 


dF 

drh 


■  (rh  -  rhn )  + 


dF 

<9x 


n 


<9x 

drh 


(rh  -  rhn) , 


(21) 


recognizing  the  dependence  of  F  upon  r/(  both  directly  as  well  as  indirectly 
through  x.  Substituting  Equations  (19)  and  (21)  into  Equation  (16),  and 
taking  into  account  the  nominal  solution  of  Equation  (17)  as  well  as  the 
approximation  of  Equation  (20), 


d  /  <9x 
dt  I  drh 


dF 

<9x 


dx 

drh 


+ 


dF 

drh 


or 


dx 


U+^ 

drh 


(22) 


Recall  that  the  vector  function  F  is  linear  in  x,  with  no  additive  constant. 
Therefore, 

clF 

“  fot  'X' 
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Substituting  into  Equation  (17), 


x 


n 


dF 


Note  the  similarity  of  this  equation  to  the  homogeneous  part  of  Equation  (22). 
This  indicates  that  any  column  of  the  homogeneous  solution  of  Equation  (22) 
takes  the  form  of  the  general  solution  of  Equation  (17). 


A. 2  Linear  Equations 

It  may  be  seen  that,  for  the  purpose  of  determining  the  gross  effects  of  hub 
position  error  on  the  relative  telescope  motion,  it  is  sufficient  to  consider 
only  the  first-order  contributions  to  the  differential  equations.  Say  that  the 
nominal  solution  for  x  is  known  in  the  form  of  a  Taylor  series.  Therefore, 
the  homogeneous  solution  to  Equation  (22)  is  also  known  in  the  form  of  a 
Taylor  series,  given  as 

e2 

Uh  —  Uho  +  tUhi  +  —  Uh2  +  •  •  •  • 

Similarly,  let  the  nominal  system  matrix  also  be  written  as  a  Taylor  series: 


Because  the  dominant  terms  of  F  are  void  of  the  components  of  r h,  the 
nonhomogeneous  part  of  Equation  (22)  is  of  higher  order,  with  Taylor  series 

<9F  _  dF  e2  9F 

drh  n~£drh  nl  +  2  Drh  n2  + 

Now,  consider  the  particular  solution  to  Equation  (22)  as  taking  the  form  of 
a  Taylor  series  as  well.  Because  the  nonhomogeneous  part  of  the  equation 
begins  at  order  1,  so  does  the  particular  solution: 

e2 

Up  —  eUpi  +  —  UP2  +  •  •  •  . 

Substituting  these  series  expansions  into  Equation  (22)  gives,  at  order 
zero, 

<9F 

Uh0  =  —  ■  Uh0.  (23) 

nO 


dF  _  <9F  <)F 

9 X  n  ~  9 X  no  +  6  9X 
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As  previously  mentioned,  this  solution  is  assumed  to  be  known  from  the 
knowledge  of  the  nominal  motion.  At  order  one, 


Uhi  = 


dF 

<9x 


•  Uhi  + 

nO 


dF 

<9x 


•  Uh  o, 

n  1 


the  solution  to  which  is  again  assumed  to  be  known,  and 


rr 

Upl  d* 


n  0 


TT  ,  9F 

•  Uv  i  +  — — 

dry. 


n  1 


(24) 


For  Equation  (24),  only  the  particular  solution  is  required.  The  system  ma¬ 
trix  here  is  the  same  as  that  of  Equation  (23);  terms  of  the  homogeneous  solu¬ 
tion  to  Equation  (24)  are  already  provided  in  the  solution  to  Equation  (23). 
It  is  assumed  that  the  solutions  are  convergent  for  any  given  time;  there¬ 
fore,  the  gross  contribution  to  the  particular  solution  may  be  found  from  this 
equation  alone.  Accordingly,  henceforth,  the  matrix  U  is  assumed  to  refer  to 
those  terms  through  Uhi  and  Uv\ . 

For  purposes  of  this  analysis,  the  order  0  terms  are  considered  to  be  those 
with  the  numerical  coefficient  A3;  order  1  terms  are  those  with  the  coefficient 
A4.  Other  terms,  including  those  involving  the  eccentricity  e  are  treated  as 
higher  order.  For  that  reason,  the  mean  motion  n  of  Earth’s  orbit  may  be 
approximated  by  nc. 


A. 3  Solution  Development 

The  solution  for  each  vector  of  the  linear  system  is  constructed  in  the  stan¬ 
dard  fashion.  Let  11  represent  a  column  vector  of  U,  and  v  represent  the  cor¬ 
responding  column  of  the  nonhomogeneous  part  of  Equation  (22),  through 
the  prescribed  order.  Then,  the  linear  system  takes  the  form 

u  =  Au  +  v, 

where  the  matrix  A  refers  to  the  matrix  of  partial  derivatives  seen  in  Equa¬ 
tion  (23).  The  homogeneous  solution  may  be  constructed  using  standard 
methods,  as  will  be  seen  below;  therefore,  it  is  necessary  to  develop  the  par¬ 
ticular  solution  to  the  system  of  differential  equations  given  by  Equation  (24). 

In  these  differential  equations,  the  nonhomogeneous  part  contains  trigono¬ 
metric  terms,  which  may  be  written  as  exponential  terms  with  exponential 
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multipliers  equal  to  the  characteristic  multipliers  of  the  homogeneous  system. 
Accordingly,  the  particular  solution  is  developed  for  the  case  of  a  single  such 
exponential  term  at  a  time;  the  solution  for  a  number  of  such  terms  may  be 
constructed  by  superposition  of  the  individual  solutions.  It  is  assumed  that 
the  characteristic  multipliers  are  distinct,  as  is  the  case  for  the  hub-telescope 
system. 

A. 3.1  Non-secular  Solution 

First,  examine  the  case  where  the  nonhomogeneous  vector  v  contains  periodic 
(or,  equivalently,  exponential)  contributions  with  frequencies  not  equal  to  one 
of  the  eigenvectors  of  the  system  matrix  A. 

Consider  a  general  first-order  nonhomogeneous  linear  system  of  differen¬ 
tial  equations  given  by 

x  =  Ax  +  elfltej, 

where  e3  is  the  jth  system  unit  vector,  and  ifl  is  not  one  of  the  system 
eigenvalues.  The  particular  solution  of  such  a  system  may  be  found  using 
the  method  of  undetermined  coefficients.  Say  that  the  particular  solution  is 
of  the  form 

xp  =  be*m, 

where  the  vector  b  is  a  constant  vector  which  is  to  be  determined.  Substi¬ 
tuting  Xp  into  the  differential  equation, 

iQ  beim  =  Aheint  +  elQter 

Equating  coefficients  of  the  exponential, 

(ifll  —  A)b  =  ej.  (25) 

As  in  is  not  an  eigenvalue  of  A,  the  solution  may  be  formed  as 

b  =  (im  —  A)~1ej. 

For  the  problem  of  interest,  the  non-secular  particular  solution  is  a  linear 
combination  of  particular  solutions  of  this  type,  each  scaled  by  the  appropri¬ 
ate  constant  coefficient  found  in  the  differential  equation.  In  the  system 

u  =  Au  +  v, 
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say  that  the  individual  contributions  to  the  particular  solution  are  now  writ¬ 
ten  as 

upi  =  db  ieinit, 

where  each  ct  is  the  appropriate  scaling  constant.  Then,  the  full  particular 
solution  is  the  sum  of  these  solutions, 


i 


Combining  with  the  homogeneous  solution,  the  complete  solution  for  u  is 

u  =  ^  (e*ni*  _  eM )  Qbj. 

i 

Note  that  this  solution  satishes  the  requirement  of  zero  initial  conditions  for 
the  variational  equations. 

A. 3. 2  Secular  Solution 

Next,  examine  the  case  where  the  nonhomogeneous  vector  v  contains  periodic 
(or  exponential)  contributions  with  frequencies  equal  to  one  of  the  eigenvalues 
of  A.  In  this  case,  a  secular  solution  will  result. 

Consider  a  general  first-order  nonhomogeneous  linear  system  of  differen¬ 
tial  equations  given  by 

x  =  Ax  +  eXitej, 

where  Xt  is  the  ith  eigenvector  of  the  matrix  A. 

The  particular  solution  of  such  a  system  may  again  be  found  using  the 
method  of  undetermined  coefficients.  Say  that  the  particular  solution  is  of 
the  form 

xp  =  (a  t  +  b)eXit, 

where  the  vectors  a  and  b  are  constant  vectors  which  are  to  be  determined. 
Substituting  xp  into  the  differential  equation, 

(a  T  A*a t  T  A;b)e  1  —  A(af  T  b)e  1  T  c  1  g j. 

Equating  coefficients  of  t  and  unity, 

(A  J  -  A)a  =  0  (26) 
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and 


(27) 


(XJ  —  .A)b  —  — a.  g^. 

Equation  (26)  may  be  solved  to  give 

a  — 

where  x,  is  the  ith  eigenvector  of  A  and  a  is  yet  to  be  determined.  Substi¬ 
tuting  into  Equation  (27), 


(A J  -  A)b  =  -ax*  +  ej. 

Next,  let  yi  denote  the  ith  left  eigenvector  of  A,  satisfying 


y?;T(Aj/  —  A)  =  0T. 


(28) 


The  elements  of  y,  describe  the  linear  dependence  of  the  rows  of  \I  —  A. 
Therefore,  a  summation  of  the  weighted  rows  of  Equation  (28)  allows  the 
elimination  of  b.  This  is  performed  by  writing 


0T  =  yiT(XiI  -  A) b  =  — y iTaxi  +  y^e,- 
=  -y  iTaxi  +  yij, 


where  y%3  is  the  jth  element  of  y,.  Solving  for  a, 


a  = 


Vij 

y  iTXj ' 


(29) 


Let  the  columns  of  A*/  —  A  be  given  by  the  column  vectors  c*.  Therefore, 
an  n-dimensional  Equation  (28)  may  be  represented  as 


[  Cl  c-2  •  •  •  cn  ]  b  =  -axi  +  ej. 

Now,  examine  all  but  the  jth  row  of  this  equation: 

[  Cl  C2  •  •  •  Cn  ]  b  =  Mb 

=  -axi, 


where  an  underbar  denotes  the  removal  of  row  j. 

The  eigenvalue  A i  is  assumed  to  be  distinct.  Therefore,  it  is  possible  to 
arbitrarily  set  the  jth  element  of  b,  bj,  to  unity,  and  solve  for  the  remainder 
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of  the  vector.  By  setting  bj  to  unity,  the  contributions  associated  with  the 
jth  column  of  \I  —  A  may  be  separated  from  the  remainder  of  the  equation 
as 

Mb  =  —  ax*  —  c  j, 

where  an  overbar  denotes  the  removal  of  column  j.  Because  the  eigenvalue  is 
distinct,  there  is  no  zero  eigenvalue,  and  a  single  row  and  column  have  been 
removed  from  \I  —  A,  the  remaining  matrix,  M,  is  nonsingular.  Therefore, 
solving  for  b, 

b  =  —  M  1  +  c j'j  .  (30) 

The  vector  b  is  then  formed  by  augmenting  b  with  the  inclusion  of  unity  as 
the  jth  element. 

To  illustrate  this  procedure,  consider  the  following  example: 


Here,  the  exponential  multiplier  associated  with  the  nonhomogeneous  part 
of  the  equation  is  1,  which  is  also  an  eigenvalue  of  the  system  matrix.  The 
corresponding  eigenvector  is 


Xi  = 


1 

1 


and  the  associated  left  eigenvector  is 

yi  =  [  3  -if. 

Equation  (29)  gives  a  as 

3  _  3 
(3)(l)  +  (-l)(l)  “2; 

Equation  (30)  gives  the  single  element  vector  b  as 


b 


-(1/3) 


|  (1)  +  (—3) 


1 

2' 


Therefore,  the  particular  solution  is 


'  1 ' 

1 

(2 

1 

t  + 

1/2 

eh 
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For  the  problem  of  interest,  the  particular  solution  is  a  linear  combination 
of  particular  solutions  of  this  type,  each  scaled  by  the  appropriate  constant 
coefficient  found  in  the  differential  equation.  In  the  system 

u  =  Au  +  v, 

say  that  the  individual  contributions  to  the  particular  solution  are  now  writ¬ 
ten  as 

upi  =  Ci  (cqtx,;  +  b,)  eXit, 

where  ct  is  the  appropriate  scaling  constant.  Then,  the  full  particular  solution 
is  the  sum  of  these  solutions, 


i 


Combining  with  the  homogeneous  solution,  the  complete  solution  for  u  is 

u  =  -eAt  ^  Qbj  +  Up. 

i 


Note  that  this  solution  satisfies  the  requirement  of  zero  initial  conditions  for 
the  variational  equations. 


A. 3. 3  State  Transition  Matrix 

In  both  the  secular  and  non-secular  cases,  the  exponential  state  transition 
matrix  is  constructed  in  the  usual  fashion.  Letting  P  denote  the  matrix  of 
eigenvectors  of  the  matrix  A,  and  A  the  diagonal  matrix  of  eigenvalues,  the 
exponential  matrix  is  given  by 


eAt  =  PAP'1. 

A. 4  Out-of-Plane  Solution 

This  approach  is  now  applied  to  the  z  and  i-components  of  Equation  (22), 
which  are  not  coupled  with  the  in-plane  components.  These  out-of-plane 
components  are 


Uz 


0  1 
-A3  0 


Uz  +  3A4 


0  0  0 

z  0  x 


(31) 
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where 


dz  dz  dz 
U  =  dxh  dyh  dzh 
dz  dz  dz 
.  dxh  dyh  dzh  J  n 

For  a  single  vector  u  of  Uz, 

u  +  v,  (32) 

where  v  is  the  appropriate  nonhomogeneous  column  of  Equation  (31). 

For  each  non-zero  column  v,  the  periodic  forcing  function  is  of  frequency 
which  is  the  same  or  approximately  the  same  as  the  natural  frequency  of  the 
system.  (If  higher-order  frequency  compensation  is  employed,  the  low-order 
frequencies  of  the  two  non- zero  columns  will,  in  fact,  be  identical.)  Therefore, 
the  secular  response  associated  with  each  non-zero  forcing  column  is  exam¬ 
ined.  Additionally,  the  non-secular  response  is  examined  for  the  case  where 
the  third-column  forcing  function  possesses  the  nearly  resonant  frequency  of 
the  linear  solution. 


u  = 


0  1 
-A3  0 


A. 4.1  Secular 


First  consider  the  secular  case,  where  v  is  of  the  form 

_  r  0 

V  Asin(y/A^t  +  ip) 

This  corresponds  to  the  first  nonhomogeneous  column  of  Equation  (31).  (If 
the  third  column  is  treated  as  secular,  it  differs  only  by  a  90°  shift  in  phase.) 
Using  exponential  representation, 


where 


v  =  A\ 


0 

A%t 


+  A2 


0 

j—iy/Ast 


A 

Ai  —  —  (sin  ip  —  i  cos  'ip) 


and 


A 

A-2  =  —(sin  ip  +  i  cos  'ip). 
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The  particular  solution  may  now  be  determined  using  the  method  of  Sec¬ 
tion  A. 3. 2,  treating  each  of  the  two  exponential  vectors  of  v  individually. 

Using  the  earlier  notation,  let  Ai  =  i\fA3  and  A2  =  —  i\[AA3.  The  eigen¬ 
vectors  are 


Xl  = 

the  left  eigenvectors  are 


1 

iVA3 


and  x2  = 


-iVA3 


yi  =  [  i  ]  and  y2  =  [  -iVA3  l  ]  . 

For  the  nonhomogeneous  contribution  given  by  e*'d^*e2,  the  method  of  un¬ 
determined  coefficients  gives 

1 

a  =  - -==, 

i2^A3 

M  =  i\/A3 , 


and 


b  = 


2A3  a/A3 
1 

Therefore,  this  part  of  the  particular  solution  is 


Upi  =  Ai 


( 

i 

■  i 

i 

\ 

1 

2 

Va3 

t  + 

2A3 

VA', 3 

V 

i 

1 

/ 

In  similar  fashion,  the  part  of  the  particular  solution  corresponding  to 
e-*vA3 1£2,  js  formed.  Combining  gives  the  full  particular  solution  as 


/ 

i 

'  1 

i 

\ 

A\ 

1 

o 

V^3 

t  + 

2A3 

\/^3 

z 

V 

1 

1 

/ 

+  A 

/ 

1 

2  c\ 

i 

VA3 

t  + 

1  i 

2A3  a /  A3 

\ 

Z 

\ 

1 

1 

/ 
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Figure  10: 


dz 

dxh 


and 


dz 

dzh 


with  resonance 


Referring  to  Section  A. 3. 2,  the  behavior  of  the  position  partial  derivative 
over  180  days  is  presented  in  Figure  10,  for  a  pair  of  fundamental  values  of  the 
phase  angle.  As  previously  mentioned,  the  hub  position  error  components  are 
taken  to  be  1  km;  therefore,  in  this  and  all  other  figures,  the  contributions  to 
the  telescope  motion  variation  is  the  same  as  the  position  partial  derivatives 
themselves. 

A. 4. 2  Non-Secular 

Next,  consider  the  response  associated  with  Equation  (32),  where  the  vector 
v  is  a  nearly-resonant  vector  of  the  form 

Acos(ut  +  0)e2. 

From  the  earlier  analysis,  the  phase  angle  0  is  taken  as  '0+90°.  The  frequency 
c o  is  taken  as  the  natural  linear  frequency  of  the  in-plane  motion.  As  shown 
in  the  earlier  analysis,  this  frequency  is  given  by 

uj  =  \Jn2  -  A3/2  +  v7 (n2  -  A3/2)2  -  (n2  +  2 A3)(n2  -  A3). 
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For  comparison,  u  is  approximately  0.035385  rad/day,  while  y/AA,  is  approx¬ 
imately  0.034148  rad/day. 

Following  the  method  of  Section  A. 3.1,  the  vector  b  is  constructed  as 

(ul  -  A)_1e2. 

The  resulting  position  partials  in  the  solution  vector  u  are  presented  in  Fig¬ 
ure  11. 

A. 5  In-Plane  Particular  Solution 

The  in-plane  ( x  and  y)  components  of  Equation  (22)  are  now  treated.  These 
differential  equations  are  of  the  form 

0  0  1  0  1  [  0  0  0  ' 

Uxy  =  2  Oo  4  n  ^  }  Uxy  +  3Aa  ^  0  0  ,  (33) 

y  n 2  +  2A3  0  0  2 n  xy  -2x  y  z  ’  v  ’ 

0  n2  —  A3  —2  n  0  y  x  0 
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where 


dx 

dx 

dx 

dxh 

dyh 

dzh 

dy 

dy 

dy 

dxh 

dyh 

dzh 

dx 

dx 

dx 

dxh 

dyh 

dzh 

dy 

dy 

dy 

dxh 

dyh 

dzh 

For  a  single  vector  u  of  Uxy, 


u  = 


0 

0 

n2  +  2A3 
0 


0 

0 

0 

n2  -  A3 


1  0 
0  1 
0  2  n 

—2  n  0 


u  +  v, 


(34) 


where  v  is  the  appropriate  column  of  the  above  matrix. 

From  the  earlier  analysis  of  the  linearized  circular  restricted  three-body 
problem,  the  system  matrix  of  Equation  (34)  has  one  eigenvalue  correspond¬ 
ing  to  a  convergent  solution,  one  corresponding  to  a  divergent  solution,  and 
a  complex  conjugate  pair  which  corresponds  to  a  periodic  solution.  In  that 
analysis,  it  was  specified  that  the  initial  conditions  would  be  selected  in  such 
a  manner  as  to  excite  only  the  periodic  modes.  Therefore,  as  in  the  previ¬ 
ous  section,  it  is  assumed  here  that  x  and  y  in  Equation  (34)  exhibit  that 
behavior;  z  has  already  been  seen  to  be  periodic. 


A. 5.1  Secular 

First,  examine  the  secular  solutions  associated  with  resonant  forcing  vectors 
v. 

Again  using  the  notation  of  Section  A. 3. 2,  consider  the  eigenvalues  given 
by  Ai  =  iu  and  A2  =  —iu>.  The  corresponding  eigenvectors  are 


i2nu> 

—i2nu 

Xl  = 

—2nu2 

and 

x2  = 

-2  nu2 

-iujx 

iujx 
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where 

X  =  <^2  T  '^2  T  2A3; 

the  left  eigenvectors  are 


yi 


i2nuj{n 2  +  2A3) 

(n2  -A3)x 

—2nuj2 

iuj\ 


and 


y2 


—i2nuj(n2  +  2A3 ) 
( n 2  ~A3)x 

—2nu2 

-iujx 


Consider  the  contribution  associated  with  the  second  column  of  Equa¬ 
tion  (33)  as  the  vector  v  of  Equation  (34).  The  order  1  solutions  for  x  and 
y  are  of  the  form 


3  A4x  =  —Acos(cjt  +  0)  and  3  A4y  =  kAsin(cjt  +  0), 


where 

k  =  ^~. 

2um 

These  functions  may  then  again  be  written  as  combinations  of  complex  ex¬ 
ponentials. 

For  the  nonhomogeneous  contribution  to  v  of  etuJte3,  the  scalar  a  is  given 


by 


a  = 


2  nuj2 

~T~ 


where 


d  =  u6  +  5(n2  +  A3)lu4  —  (5  n2  —  4  A3)(n2  +  2A3)lu2  —  (n2  —  A3)(n2  +  2  A3)2. 


The  associated  vector  b  is 


4  n2u2/d  —  i/uj 

2n/C  +  AnxuJ3  /  (d() 

1 

—2  nu2x(u2  —  n2  +  A3)/(d()  +  i2nu/C, 
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Figure  12: 


dx  dy 
and 


d'Uh  dyh 


with  resonance 


where 

C  =  U2  +  n2  -  A3. 

For  the  contribution  of  e~luJte3,  the  scalar  a  is  the  same,  and  the  associated 
vector  b  is  the  complex  conjugate. 

For  the  contribution  of  elute4, 


and 

—2 n/x  ~  i^nuj3 /d 

X2/d  -  i/u 

2nuj2{uj2 —  n2 —  2A3)/d  —  i2nuj/x 

1 

For  the  contribution  of  elujte^  the  resulting  a  and  b  are  the  complex  conju¬ 
gates  of  these. 

Combining  into  the  solution  vector  u,  the  resulting  position  partial  deriva¬ 
tives,  in  the  second  column  of  matrix  Uxy,  are  presented  in  Figures  12(a) 
(^-component)  and  12(b)  (^/-component). 
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Figure  13: 


dx  ,  dy 
and 


dx 


dxh 


with  resonance 


Next,  the  same  procedure  is  followed  for  the  first  column.  The  solution  is 
constructed  using  the  same  a  scalars  and  b  vectors  as  above.  The  resulting 
position  partial  derivatives  are  presented  in  Figures  13(a)  (x- component) 
and  13(b)  (^/-component). 

Finally,  the  procedure  is  followed  for  the  third  column.  If  the  low-order 
solution  for  z  is  taken  as  being  resonant  with  the  system  matrix,  the  solution 
associated  with  the  forcing  function  SA^ze^  may  be  constructed  using  the 
same  a  and  b  expressions  as  above.  The  resulting  position  partial  derivatives 
are  presented  in  Figures  14(a)  (x-component)  and  14(b)  (^/-component). 

A. 5. 2  Non-Secular 

Finally,  the  case  is  treated  where  the  uncompensated  low-order  solution  for  z 
is  considered  to  be  non-resonant  with  the  in-plane  system.  The  third-column 
forcing  vector  is  taken  to  be  of  the  form 

Ask^v^M  +  0)e3- 

Again  following  the  method  of  Section  A. 3.1,  the  vector  b  is  constructed  as 

(oT  -  A)_1e3. 
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The  resulting  solution  vector  u  is  presented  in  Figures  15(a)  (x- component) 
and  15(b)  (^-component). 

A.  6  Summary 

From  the  plots  of  this  appendix,  it  is  recognized  that  all  of  the  small- 
magnitude  contributions  to  the  solution  of  the  variational  equations  are  in 
the  millimeter  range.  Over  180  days,  the  large-magnitude  contributions  reach 
hundreds  of  millimeters.  However,  even  these  contributions  remain  extremely 
small  over  the  first  90  days  from  epoch.  Therefore,  it  is  concluded  that  the 
errors  which  are  induced  by  1  km  errors  in  hub  position  are  sufficiently  small 
that  they  may  be  ignored. 
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B  Relative  Motion  Near  Earth-Moon  L2  Li- 
bration  Point 

At  the  sponsor’s  request,  a  rudimentary  analysis  was  performed  into  the 
approach  which  would  be  taken  in  an  investigation  of  relative  motion  near 
the  Earth-Moon  L2  point.  In  this  system,  Earth  and  Moon  are  the  pri¬ 
mary  bodies;  here,  the  fundamental  periods  in  and  out  of  the  xy-plane  are 
approximately  14.650  days  and  15.275  days,  respectively. 

Consider  the  case  of  a  hub  orbiting  the  L2  point  at  roughly  a  distance 
of  30,000  km,  with  a  telescope  located  100  m  from  the  hub.  The  effects  of 
various  terms  of  the  elliptical  restricted  expansion  over  15  days  are  given  in 
Table  5. 

Table  5:  Along-Ellipsoid  Effect  of  Earth-Moon  Perturbation  on  Solution  (15 
days) 


perturbing 

term 

effect 

(m) 

re 

21.5 

rrh 

114.1 

rrhe 

8.4 

rrh 2 

64.4 

rrh 3 

31.5 

4 

rrh 

14.9 

The  effects  are  shown  for  terms  which  contribute  to  as  low  as  approxi¬ 
mately  10  m  over  the  15  days.  It  is  noted  that  the  largest  contributions  from 
perturbations  by  Sun  during  this  time  period  are  on  the  order  of  less  than 
1  m.  Perturbing  forces  from  solar  radiation  pressure  and  spacecraft  thruster 
firing  are  on  the  same  order  as  for  the  Sun-Earth  system. 


164 


C  Calculation  of  Luminosity  Reduction  Due 
to  Partial  Eclipse 

The  calculation  of  the  luminosity  reduction  factor,  er,  is  shown  here.  This 
symbol  is  a  variable  within  Equation  (8),  which  is  used  to  compute  the  force 
per  unit  mass  due  to  solar  radiation  pressure. 


Figure  16:  Solar  Radiation  Shadowing  Geometry 

Figure  16  shows  the  shadowing  geometry.  In  the  figure,  line  CS  is  normal 
to  line  PQ;  line  CD  is  normal  to  the  Sun-Earth  line,  SE;  line  SO  is  normal 
to  line  CO.  The  distance  Rj,  is  the  Earth’s  effective  obscuration  radius.  Dis¬ 
tances  Re  and  Rs  are  the  mean  radii  of  Earth  and  Sun,  respectively.  Solar 
shadowing  is  based  on  projecting  2/4  onto  line  PQ. 

The  calculation  of  a  is  a  matter  of  the  geometry  of  two  intersecting  circles. 
Figure  17  shows  the  additional  geometry  and  the  notation  used  in  the  cal¬ 
culation.  The  projection  of  the  obscuration  disk  on  the  plane  normal  to  the 
spacecraft-Sun  line  (CS)  is  an  ellipse  of  negligible  eccentricity.  This  ellipse 
is  very  closely  approximated  as  a  circle.  The  comment  lines  in  the  accom¬ 
panying  FORTRAN  code  describes  the  notations  used  in  Figures  16  and  17. 
The  routine  calculates  the  solar  radiation  force  per  unit  mass  at  any  distance 
greater  than  the  Sun-Earth  distance.  The  code  also  calculates  the  diminished 
radiation  force  for  all  spacecraft  trajectories  that  pass  through  any  portion 
of  the  Earth’s  shadow.  It  is  assumed  that  the  radiation  force  falls  off  linearly 
in  proportion  to  the  percentage  of  solar  area  lost  to  obscuration. 
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Figure  17:  Earth  Obscuration  and  Solar  Disk  Geometry 


PROGRAM  SOLRAD 
C 

IMPLICIT  REAL*8 (A-H , 0-Z) 

REAL*8  LRF , MASS 
COMMON/INPUTS/CR , AREA , MASS 
C 

C  EXAMPLE  DRIVER  FOR  SOLAR  RADIATION  FORCE 

C 

CR=1 . 5D0 
AREA=10 . DO 
MASS=2000 . DO 
SE=149597870 . DO 
C  SEE  FIG  16: 

WRITE (*,*)  J  INPUT  CE  AND  CD  =  J 
READ ( * , * )  CE,CD 
C 

CALL  LRF  ACTOR ( SE , CE , CD , LRF ) 

WRITE (*,*)  J  LRF  =  ' , LRF 
C 

CS=DSQRT(SE**2+CE**2) 

CALL  RFORCE (LRF , CS , F) 

WRITE(* , 101)  F 

101  FORMAT ( ’  SOLAR  RADIATION  FORCE/MASS  (M/S~2)  =  ’ .1PD11.4) 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


STOP 

END 

SUBROUTINE  LRFACTOR (SE , CE , CD , LRF) 

IMPLICIT  REAL*8 (A-H , O-Z) 

REAL*8  LRF,L 
DATA  RS/ 696000 . DO/ 

DATA  RE/6378. 14D0/ 

DATA  PI/3. 141592653589793D0/ 

LRF  =  LUMINOSITY  REDUCTION  FACTOR  (OUTPUT) 

FIG  16  CALCULATIONS: 

RS  =  SOLAR  RADIUS  (KM) 

RE  =  EARTH  RADIUS  (KM) 

SE  =  SUN  EARTH  DISTANCE  (KM)  (INPUT) 

CE  =  SATELLITE  TO  EARTH  DISTANCE  (KM)  (INPUT) 

CE  IS  A  PORTION  OF  CO. 

CD  =  NORMAL  DISTANCE  OF  SATELLITE  TO  SUN-EARTH 
LINE  (KM)  (INPUT)  (GE  ZERO) 

CS  =  SATELLITE  TO  SUN  DISTANCE  (KM) 

CO  =  SATELLITE  TO  SUN  (THRU  E)  OFFSET  DISTANCE  (KM) 

OS  =  PERPENDICULAR  DISTANCE  OF  SUN  FROM  CO  (KM) 

CO  +  OS  +  CS  FORMS  A  RIGHT  TRIANGLE  WITH 
HYPOTENUSE  CS . 

PS  =  NORTHERN  LIMIT  OF  SOLAR  OSBSCURATION  ALONG  LINE  PQ 
NORMAL  TO  CS .  PS  MEASURED  FROM  SOLAR  CENTER  S. 

FIG  17  CALCULATIONS: 

RO  =  RADIUS  OF  OBSCURATION  DISK  (KM) 

C  =  SEPARATION  BETWEEN  GEOMETRIC  CENTER  OF  SOLAR  DISK 
AND  GEOMETRIC  CENTER  OF  OBSCURATION  DISK  (KM) 

=  SMALLA  +  SMALLB 
=  SO  ’ 
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C  L  =  DISTANCE  FROM  LINE  SEGMENT  C  TO  EITHER  POINT  OF 

C  INTERSECTION  OF  THE  2  DISK  PERIMETERS  (KM) 

C  L**2  =  R0**2  -  SMALLA**2  =  Rl**2  -  SMALLB**2  (RIGHT 

C  TRIANGLES) 

C  THETAO  =  DATAN (L/SMALLA) 

C  THETA 1  =  DATAN (L/SMALLB) 

C  SO  =  SECTOR  AREA  SUBTENDED  BY  2*THETA0  (KM**2) 

C  SI  =  SECTOR  AREA  SUBTENDED  BY  2*THETA1  (KM**2) 

C  ASHADOW  =  AREA  OF  INTERSECTION  OF  OBSCURATION  DISK  AND 

C  SOLAR  DISK  (KM**2) 

C  SUNAREA  =  AREA  OF  SOLAR  DISK  (KM**2) 

C 

C  SEE  FIG  16  FOR  THE  FOLLOWING: 

C 

ED=DSQRT (CE**2-CD**2) 

C  COSINE (LAMBDA) 

COSL=ED/CE 
CO=SE*COSL+CE 
C  SINE (LAMBDA) 

SINL=CD/CE 

OS=SE*SINL 

CA=DSQRT (CE**2-RE**2) 

CS=DSQRT  ( (SE+ED)  **2+CD**2) 

C  OBSCURATION  OCCURS  WHEN  DABS (PS) . LE . RS 

PS=CS* (RE*CO-OS*CA) / (CA*CO+OS*RE) 

IF (PS . LT . -RS)  THEN 
C  NO  ECLIPSE 

LRF=1 . ODO 
RETURN 
END  IF 

IF (PS . GT . RS)  THEN 
C  TOTAL  ECLIPSE 

LRF=0 . ODO 
RETURN 
END  IF 
C 

C  PARTIAL  ECLIPSE  (DABS (PS) .LE.RS) 

C 
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C  COSINE (BETA  -  ALPHA) 

COSBA=CO/CS 

C  "NORTHERN"  RADIUS  OF  OBSCURATION  MEASURED  FROM  O’ 

C  "NORTH"  TO  P  ALONG  LINE  PQ : 

SOP=OS/COSBA 

RADIUSN=PS+SOP 

RO=RADIUSN 

C 

C  SEE  FIG  17  FOR  THE  FOLLOWING: 

C 

C=SOP 

SMALLB= (RS**2-R0**2+C**2) / (2 . DO*C) 

SMALLA=C-SMALLB 

L=DSQRT(RS**2-SMALLB**2) 

THETAO=DATAN (L/SMALLA) 

THETA1=DATAN (L/SMALLB) 

S0=R0**2*THETA0 
S1=RS**2*THETA1 
ASHADOW=SO+Sl-C*L 
SUNAREA=PI*RS**2 
LRF=1 . DO-ASHADOW/SUNAREA 
RETURN 
END 
C 

SUBROUTINE  RFORCE (LRF , CS , F) 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  LRF, MASS 
COMMON/INPUTS/CR , AREA , MASS 
C 

C  LRF  =  LUMINOSITY  REDUCTION  FACTOR  (INPUT) 

C  CS  =  HELIOCENTRIC  DISTANCE  TO  NGST  (KM)  (INPUT) 

C  F  =  SOLAR  RADIATION  FORCE  PER  UNIT  MASS 

C  (NEWTONS/KG  =  M/S~2)  (OUTPUT) 

C 

C  CR  =  SOLAR  FLUX  REFLECTION  PARAMETER  0  <=  CR  <=  2 

C  (INPUT  THRU  COMMON) 

C  AREA  =  NGST  AREA  PROJECTED  NORMAL  TO  SUN  LINE  CE 

C  (METERS* *2) 
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C  (INPUT  THRU  COMMON) 

C  MASS  =  NGST  MASS  (KG)  (INPUT  THRU  COMMON) 

C 

F=1 . 0198D17*CR*AREA*LRF/ (1000 . DO*CS) **2/MASS 

RETURN 

END 
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D  Constant  Parameters 


Table  6:  Physical  Constants  [5] 


gravitational  parameter, 

Sun  alone 

AC 

132,712,440,017.987  km3/s2 

gravitational  parameter, 
Earth  alone 

^2 

398,600.4415  km3/s2 

gravitational  parameter, 
Earth-Moon 

^2 

403,503.236  km3/s2 

gravitational  parameter, 
Moon  alone 

AU 

4,902.8003  knr3/s2 

astronomical  unit 

AU 

149,597,870.691  km 

mean  Earth-Moon  barycen- 
ter  distance  from  Sun 

1.000001018  AU 

eccentricity  of  Earth-Moon 
barycenter  orbit  about  Sun 

e 

0.01670862 

mean  motion  of  Earth-Moon 
barycenter  orbit  about  Sun 

n 

0.199106385  xl0“6  rad/s 

L2  distance  ratio 

7 

0.01007824 

Table  7:  Computed  Values  and  Coefficients 


reference  Sun- 
Earth  distance 

D 

149,618,905.218739  km 

gravitational 

A3 

1.16556055765939  x  10“3  1/day2 

coefficients 

a4 

5.84525170441422  x  10~10  l/(km-day2) 

(see  page  126) 

A5 

3.86396147244215  x  10~16  l/(km2-day2) 

A6 

2.56240418152728  x  10"22  l/(km3-day2) 
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